Historical Demographic Processes Dominate Genetic Variation in Ancient Atlantic Cod Mitogenomes

Martínez-García, Lourdes; Ferrari, Giada; Oosting, Tom; Ballantyne, Rachel; van der Jagt, Inge; Ystgaard, Ingrid; Harland, Jennifer; Nicholson, Rebecca; Hamilton-Dyer, Sheila; Baalsrud, Helle Tessand; Brieuc, Marine Servane Ono; Atmore, Lane M.; Burns, Finlay; Schmölcke, Ulrich; Jakobsen, Kjetill S.; Jentoft, Sissel; Orton, David; Hufthammer, Anne Karin; Barrett, James H.; Star, Bastiaan

doi:10.3389/fevo.2021.671281

ORIGINAL RESEARCH article

Front. Ecol. Evol., 04 June 2021

Sec. Evolutionary and Population Genetics

Volume 9 - 2021 | https://doi.org/10.3389/fevo.2021.671281

This article is part of the Research Topic Assessing Biodiversity in the Phylogenomic Era View all 10 articles

Historical Demographic Processes Dominate Genetic Variation in Ancient Atlantic Cod Mitogenomes

$\r\nLourdes Martínez-García*&#x;$ Lourdes Martínez-García^1*†

Giada Ferrari^1†

Tom Oosting^2†

Rachel Ballantyne^3†

Inge van der Jagt⁴

Ingrid Ystgaard^5,6

Jennifer Harland⁷

Rebecca Nicholson⁸

Sheila Hamilton-Dyer⁹

Helle Tessand Baalsrud^1†

Marine Servane Ono Brieuc^1†

Lane M. Atmore^1†

Finlay Burns¹⁰

Ulrich Schmölcke¹¹

Kjetill S. Jakobsen^1†

Sissel Jentoft^1†

David Orton^12†

Anne Karin Hufthammer¹³

James H. Barrett^3,6†

Bastiaan Star^1*†

¹Department of Biosciences, Centre for Ecological and Evolutionary Synthesis (CEES), University of Oslo, Oslo, Norway
²School of Biological Sciences, Victoria University of Wellington, Wellington, New Zealand
³McDonald Institute for Archaeological Research, Department of Archaeology, University of Cambridge, Cambridge, United Kingdom
⁴Department of Archaeology, Cultural Heritage Agency of the Netherlands, Amersfoort, Netherlands
⁵Department of Historical and Classical Studies, Norwegian University of Science and Technology (NTNU), Trondheim, Norway
⁶Department of Archaeology and Cultural History, Norwegian University of Science and Technology (NTNU) University Museum, Norwegian University of Science and Technology, Trondheim, Norway
⁷Archaeology Institute, University of the Highlands and Islands, Orkney, United Kingdom
⁸Oxford Archaeology, Oxford, United Kingdom
⁹Department of Archaeology & Anthropology, Bournemouth University, Poole, United Kingdom
¹⁰Marine Scotland Science, Aberdeen, United Kingdom
¹¹Centre for Baltic and Scandinavian Archaeology (ZBSA), Schleswig, Germany
¹²BioArCh, Department of Archaeology, University of York, York, United Kingdom
¹³Department of Natural History, The University Museum, University of Bergen, Bergen, Norway

Ancient DNA (aDNA) approaches have been successfully used to infer the long-term impacts of climate change, domestication, and human exploitation in a range of terrestrial species. Nonetheless, studies investigating such impacts using aDNA in marine species are rare. Atlantic cod (Gadus morhua), is an economically important species that has experienced dramatic census population declines during the last century. Here, we investigated 48 ancient mitogenomes from historical specimens obtained from a range of archeological excavations in northern Europe dated up to 6,500 BCE. We compare these mitogenomes to those of 496 modern conspecifics sampled across the North Atlantic Ocean and adjacent seas. Our results confirm earlier observations of high levels of mitogenomic variation and a lack of mutation-drift equilibrium—suggestive of population expansion. Furthermore, our temporal comparison yields no evidence of measurable mitogenomic changes through time. Instead, our results indicate that mitogenomic variation in Atlantic cod reflects past demographic processes driven by major historical events (such as oscillations in sea level) and subsequent gene flow rather than contemporary fluctuations in stock abundance. Our results indicate that historical and contemporaneous anthropogenic pressures such as commercial fisheries have had little impact on mitogenomic diversity in a wide-spread marine species with high gene flow such as Atlantic cod. These observations do not contradict evidence that overfishing has had negative consequences for the abundance of Atlantic cod and the importance of genetic variation in implementing conservation strategies. Instead, these observations imply that any measures toward the demographic recovery of Atlantic cod in the eastern Atlantic, will not be constrained by recent loss of historical mitogenomic variation.

Introduction

Continuous human activities and a changing climate have influenced terrestrial and marine ecosystems for millennia (Venter et al., 2016; Rodrigues et al., 2019; Mitchell and Rawlence, 2021), impacting the evolutionary potential and population demography of a range of species (Seersholm et al., 2018). Ancient mitochondrial DNA (mtDNA) has been widely used to understand long-term genomic consequences of such impacts (Shapiro et al., 2004; Nyström et al., 2006; Stiller et al., 2010; Paijmans et al., 2013; Fortes and Paijmans, 2015; Casas-Marce et al., 2017). Nonetheless, most ancient mtDNA studies have focused on terrestrial species, and studies that investigate the impacts of long-term human activities and/or climatic variation on fish, using whole genome sequencing approaches, are relatively rare. Long-term commercial fisheries—covering many centuries—have contributed to the decline of economically and ecologically important marine species (Exadactylos et al., 2007; Pinnegar and Engelhard, 2008; Barrett, 2019). The consequences of intensive fishing in recent times may be difficult to assess as this requires an understanding of historical population dynamics (Selim et al., 2016). The analysis of long-term biological and demographic fluctuations can therefore help to improve guidelines for sustainable fisheries management and optimal conservation measures (Barrett, 2019). In order to provide a long-term perspective on fishing exploitation impacts, the use of archeological evidence, such as fish bone remains, is essential for those periods for which little or no historical data are available. Recent developments in whole genome aDNA methods now allow the inference of demographic histories and the estimation of genetic fluctuations over time from fishbone samples (Oosting et al., 2019; Ferrari et al., 2021). Such combined molecular analyses of historical and modern samples can potentially provide an understanding of the association between human-environmental impact and population declines (Hofman et al., 2015).

Several studies have shown the utility of temporal mtDNA analyses in the marine environment. For instance, ancient mitogenomes have investigated impacts of climate and hunting on the Atlantic walrus (Star et al., 2018; Keighley et al., 2019; Barrett et al., 2020), narwhals (Louis et al., 2020), and the extinct great auk (Thomas et al., 2019). In fish, such studies remain limited to partial mitogenome data. For example, a shift in sturgeon species distributions was detected during the Holocene in the North East Atlantic based on CytB amplicon data (Nikulina and Schmölcke, 2016). Moreover, impacts of habitat destruction and human activities during the 1800s were associated with a reduction of the mtDNA diversity of Chinook salmon from the Columbia River in the 12S and control region by comparing ancient and modern samples (Johnson et al., 2018). Similarly, impacts of human exploitation and climate oscillations were associated with losses of haplotypic CytB variation in Atlantic cod during the 15th to 16th centuries in Iceland (Olafsdottir et al., 2014). In contrast, comparable levels of ancient mtDNA genetic diversity were found between ancient and modern samples of herring specimens, despite continuous human exploitation (Speller et al., 2012). Notwithstanding these examples, human-environmental impacts and population declines remain unclear for a wide range of marine species and populations.

Atlantic cod (Gadus morhua L. 1758) is a benthopelagic predatory fish with high reproductive rates and with a fundamental ecological role in marine ecosystems (Barth et al., 2017; Edvardsson et al., 2019). It has been one of the most exploited fish species in the North Atlantic Ocean (Carr et al., 1995; Árnason et al., 2000; Nicholls et al., 2021). The distribution of this species extends through the cold waters of North America, across the continental shelves of Greenland and Iceland, and northern Europe (Lait et al., 2018). Relatively large population sizes have been characteristic throughout its entire distribution even during the expansion of long-distance fish trading during the 12th to 13th centuries in the eastern Atlantic and at the beginning of the 16th century in the western Atlantic (Barrett et al., 2004, 2011; Orton et al., 2014; Castañeda et al., 2020). However, intensive fishing activities during the 20th century (Mieszkowska et al., 2009; Jonsson et al., 2016; Brattey et al., 2018) resulted in the severe depletion of several stocks, for instance the North Sea stock, which was decimated from annual landings of 354,000 to 50,000 tons during this period (Bannister, 2004). In addition to past human exploitation, climatic events like the Little Ice Age—a cooling period that varied regionally in timing and duration but occurred between ca. 1300–1850 CE—may have caused large declines between the sixteenth and 17th centuries (Edvardsson et al., 2019).

The genomic consequences of such population dynamics and declines in Atlantic cod remain unclear. Based on partial and whole mtDNA data, Atlantic cod populations between the western and eastern Atlantic Ocean show significant structure (Árnason, 2004; Jørgensen et al., 2018; Lait et al., 2018), whereas low to no mtDNA differentiation has been found across a wide range of eastern Atlantic locations (Carr et al., 1995; Árnason and Palsson, 1996; Árnason et al., 1998, 2000; Sigurgíslason and Árnason, 2003). Here, we compared modern and ancient Atlantic cod mitogenomes—dated up to 6500 BCE—from different fishing locations in northern Europe. We evaluated whether Atlantic cod in the eastern Atlantic have experienced any loss of genetic variation, analyzed long term patterns of effective population size, and related any observed decline to the impact of commercial fisheries or climate change.

Materials and Methods

Sample Collection

Ancient samples of Atlantic cod (n = 48) were obtained from 11 excavation sites (Figure 1 and Supplementary Table 1) and were stored dry and unfrozen. The specimens were all supplied by the relevant archeological organizations, or sampled with permission on their premises. The shipment of Atlantic cod bones does not require CITES or other wildlife regulation permits for transport or analysis. Where practicable, only a subsample of bone was employed for the aDNA research, leaving material for other studies. Dating of the samples (Supplementary Table 1) was based on archeological context. Ancient samples were morphologically and genetically identified as Atlantic cod. A total of 472 available modern mitogenomes were obtained from Jørgensen et al. (2018), Lait et al. (2018), and Barth et al. (2019). Novel mtDNA sequence data from modern specimens sampled in 2016 in Orkney, United Kingdom (n = 24) were also included (Figure 1 and Supplementary Table 2). The collection of the Orkney specimens complied with the Nagoya Protocol and Convention on Biological Diversity, which the United Kingdom signed up to in 2016. All specimens were deceased when the fin clip was collected.

FIGURE 1

Figure 1. (A) Ancient and modern Atlantic cod sampling locations. The abbreviated name of archeological sites is shown in the map (see also Supplementary Table 1). Modern data were obtained from Jørgensen et al. (2018); Lait et al. (2018) and Barth et al. (2019). Sample sites are grouped according to larger oceanographic location, ecotype and modern or ancient specimens: western Atlantic (WA) = dark green, Baffin Island (BI) = light green, Irish Sea (IS) = gray, North Sea (NS) = purple, western Baltic (WB) = orange, eastern Baltic (EB) = yellow, Tvedestrand fjord (TV) = pink, North East Arctic (NE) = dark blue, Norwegian coast (NC) = blue, and Ancient specimens (AN) = red (i.e., Saevarhelleren = SAE, Ørland Kampflybase = ØRL, Bjørkum = BJØ, Haithabu = HAI, Quoygrew = QUO, Schleswig Schild = SCH, Oslo Mindets Tomt = OSL, Bristol Dundas Wharf = DUW, Skonsvika = SKO, King’s Lynn Raynham House = KIN, Kongshavn = KON, London Trig Lane = TRI, Vlaardingen Gat in de Markt = VLA, Bristol Finzel’s Reach = FIN, Rotterdam, Hoogstraat 13–26 = ROT, Aberdeen = ABE, Cambridge Grand Arcade = CAM, Skriðuklaustur = SKI, Newport Ship = NEW). Modern distribution of Atlantic cod is shown in light blue. Fish illustration drawn by Geir Holm. (B) Cleveland plot showing the age range of ancient samples. Sample age is estimated based on archeological context. The number of individuals per site is represented in brackets after the site name (Supplementary Table 1).

DNA Extraction, Amplification and Sequencing

DNA extraction and library preparation from ancient samples were performed in the aDNA laboratory at the University of Oslo under rigorous measures (Cooper and Poinar, 2000; Gilbert et al., 2005). All ancient samples were processed with the same DNA extraction and library protocols according to Ferrari et al. (2021). In short, bones were UV-treated for 10 min per side and pulverized using a stainless-steel mortar (Gondek et al., 2018). Per specimen, two aliquots containing between 150 and 200 mg of bone powder were used as starting material for DNA extraction. Double-indexed blunt-end sequencing libraries were built from 15 to 16 μl of DNA extract using the Meyer-Kircher protocol (Meyer and Kircher, 2010; Kircher et al., 2012) with the modifications listed in Schroeder et al. (2015) and the single-tube (BEST) protocol (Carøe et al., 2018) with the modifications described in Mak et al. (2017). Sequencing reads were processed using PALEOMIX v1.2.13 (Schubert et al., 2014). Trimming of residual adapter contamination, filtering and collapse of reads was done using AdapterRemoval v.2.1.7 (Lindgreen, 2012). Sequencing reads shorter than 25 bp were discarded. Mapping of remaining reads was performed against the Atlantic cod GadMor3.0 nuclear genome (RefSeq assembly accession GCF_902167405.1; Star et al., 2011; Tørresen et al., 2017) and mitochondrial genome (Johansen and Bakke, 1996) using BWA v.0.7.12 (Li and Durbin, 2009) with the aln algorithm, disabled seeding and minimum quality score of 25. The resulting BAM files were indexed with samtools v.1.9 (Li et al., 2009) and DNA postmortem damage assessed using MapDamage v.2.0.9 (Jónsson et al., 2013). DNA from modern Orkney samples were extracted using a DNeasy Blood & Tissue kit (Qiagen). Libraries were assembled with a TrueSeq DNA PCR-Free Preparation Kit and sequenced on an Illumina HiSeq 2,500. Modern alignment—including Orkney and Barth et al. (2019) samples, and the outgroup Alaska pollock (Gadus chalcogrammus; Malmstrøm et al., 2016) – was carried out using BWA v.0.7.12 with the mem algorithm, and a minimum quality score of 25.

Mitogenomic Analysis

Variant calling was performed using GATK v.4.1.4. (McKenna et al., 2010) simultaneously in all ancient, modern Orkney and Barth et al. (2019) samples, including the outgroup. gVCF files were created for each sample using GATK HaplotypeCaller (with ploidy set to 1). Individual genotypes were combined in one file using GATK CombineGVCFs and GenotypeGVCFs. Filtering was performed using bcftools v.1.9. (Li et al., 2009) and vcftools v.0.1.16. (Danecek et al., 2011) with the following thresholds: FS < 60.0, SOR < 4, MQ > 30.0, QD > 2.0, SnpGap = 10, minGQ = 15, minDP = 3, remove indels = yes, biallelic loci = yes, meanDP < 30 and read depth > 3. Consensus sequences were built using bcftools consensus and aligned using MAFFT v.7.429 (Katoh and Standley, 2013). Available modern mitogenomes obtained from Jørgensen et al. (2018) and Lait et al. (2018) were manually inspected using MEGA v.7 (Kumar et al., 2016) to set as missing the consistent nucleotide differences (between 50 and 100%; Supplementary Table 3) between their Illumina, Sanger and/or Roche 454 technologies with the Illumina sequenced mitogenomes in this study. Control region and half of the tRNA-Pro sequences from all the mitogenomes were excluded from further analyses as these two regions were not fully complete (i.e., 15,696–15,815 positions) after aligning sequences obtained from Jørgensen et al. (2018) and Lait et al. (2018) with the sequences presented in this study and Barth et al. (2019) samples. Thus, all sequences analyzed had 15,695 bp in length. Validated SNPs were annotated as transversion and/or transition using SNP-sites (Page et al., 2016). Checked and modified modern sequences (Jørgensen et al., 2018; Lait et al., 2018) were added and aligned to our multi-fasta alignment using MAFFT v.7.429. Unique sequences were identified with IQTREE v.1.6.12 (Nguyen et al., 2015). File formats required for different software and/or packages were obtained with seqinr and ape (i.e., nexus format; Paradis and Schliep, 2019; Charif et al., 2020), and phylotools (i.e., phylip format; Zhang et al., 2017) packages implemented in R.

Different sample combinations were used to compare the genetic diversity of the ancient samples to those of the modern conspecifics. Given the low spatial structure in the eastern Atlantic region (Árnason and Palsson, 1996; Árnason et al., 1998; Sigurgíslason and Árnason, 2003) and lack of consistent spatial structure amongst specimens (Supplementary Figures 3, 5, 6), all 48 ancient samples were compared as a single group to modern samples grouped into larger marine locations (according to their geographical proximity or ecotype; Figure 1). In addition, a comparison of subsets of multiple specimens from two archeological locations (Quoygrew and Haithabu) for which a more specific temporal pair from the same geographical region could be identified, was performed (Supplementary Table 1). Quoygrew specimens were locally sourced (Harland and Barrett, 2012; Star et al., 2017). Therefore, modern specimens sampled in the same area (i.e., modern Orkney) provide a logical, spatially consistent temporal comparison. However, specimens from Haithabu, were sourced from northern Norway (Star et al., 2017), and belonged to the North East Arctic ecotype. For these traded specimens, the North East Arctic ecotypes provide a spatially relevant temporal comparison, rather than North Sea or western Baltic specimens.

Haplotype (h) and nucleotide diversities (π), number of haplotypes (Nh) and number of polymorphic sites (S) were calculated using DnaSP v.6 (Rozas et al., 2017). To allow direct comparison with earlier CytB results (Árnason, 2004; Olafsdottir et al., 2014; Jørgensen et al., 2018), specific CytB haplotypes based on 250 bp gene fragment as previously reported by Árnason (2004) were identified using MEGA v.7. Demographic histories were determined by Tajima’s D (TD) and Fu’s F (F) neutrality in DnaSP v.6. A different number of specimens were obtained for ancient and modern locations. We corrected for such differences in sample size by randomly downsampling the modern specimens for each of the temporal comparisons (North East Arctic and Orkney) using 1,000 bootstrap replicates. A 95% confidence interval of the genetic parameters; genetic variation (π) and patterns of population demography (TD and F) was calculated from these 1,000 bootstrap replicates that were sampled using a without replacement approach with the sample function implemented in R (R Core Team, 2020) and the fasta.sample function in the FastaUtils package also in R (Salazar, 2020). For the bootstrapping test, π, TD and F from temporally spaced modern locations were re-calculated with the pegas (Paradis, 2010) and PopGenome (Pfeifer et al., 2020) packages implemented in R. Relationships among ancient and modern samples were visualized for whole mitogenome and CytB sequence data, by constructing a mitochondrial haplotype-genealogy graph using Fitchi (Matschiner, 2016) with the ML-based phylogenetic tree obtained with IQTREE v.1.6.12 as input.

Population Dynamics and Demographic Reconstruction

An analysis of molecular variance (AMOVA, 1,000,000 permutations) and population pairwise genetic distances (Φ_ST) were obtained in Arlequin v.3.5 (Excoffier and Lischer, 2010), to determine the distribution of variation between marine locations and temporally spaced locations. Divergence and coalescent analyses were based on unique sequences only (n = 525 sequences including the outgroup). Substitution model selection for unique sequences was performed using PHYML v.3.1 (Guindon et al., 2010) as implemented in JMODELTEST v.2.1.10 (Guindon and Gascuel, 2003; Darriba et al., 2012). Model selection was determined on the following partitions: 1st, 2nd, and 3rd codons from protein coding regions, rRNAs and tRNAs. Best-fitting models were selected according to the Aikaike Information Criterion (AIC; Supplementary Table 4). Based on these results, phylogenetic estimates were obtained using BEAST v.2.6.3 (Bouckaert et al., 2019).

Bayesian settings for all phylogenetic analyses included two sets of partitions: coding region and non-coding region. Three independent runs to test for chain convergence were run under the Coalescent Constant Population Tree Prior. Tip ages (ancient and modern dates) were included for each set of runs (Supplementary Tables 1, 2). Sample dates for ancient specimens were rounded to a midpoint date—from a given range—where necessary. To achieve high effective sample sizes (ESS = ≥ 200), chain lengths were run 800,000,000 under a substitution rate of 1.14 × 10^–8 substitution/site/year as per Lait (2016) assuming a GTR + I (for coding regions) and TIM1 + I (for non-coding regions) models of evolution and a strict clock. Tracer v.1.71 (Rambaut et al., 2018) was used to check for convergence of MCMC and to ensure sufficient sampling. Consensus trees were obtained using TreeAnnotator v.2.6.2—implemented in BEAST v.2.6.3—after 10% burn-in. Final phylogenetic trees were viewed and edited in FigTree v.1.4.4.

Finally, a Coalescent Bayesian Skyline (CBS) analysis was completed to reconstruct the demographic history—including female effective population size (Ne) —of Atlantic cod through time. To assess any confounding effect of past or contemporary population structure (Heller et al., 2013), we analyzed demographic history using 6 different data sets (excluding the outgroup): (I) all 524 sequences, (II) 476 modern sequences (excluding 48 ancient samples), (III) 273 sequences (excluding clades associated with most western Atlantic and Baltic Sea samples), (IV) 368 sequences (excluding the clade associated with most Baltic Sea samples), (V) 429 sequences (excluding clades associated with most western Atlantic samples) and (VI) 48 ancient sequences (excluding all modern samples). The specific clades that were excluded in III, IV and V can be found in Supplementary Figure 4. We used the same MCMC sampling procedure described before with 3 independent runs reaching convergence at high effective sample sizes (ESS = ≥ 200). Chain lengths were run 800,000,000 for data sets I, II and V with a number of bPopSize and bGroupSize of 10; while chain length for data sets III and IV were run 500,000,000 and 50,000 for data set VI with a number of bPopSize and bGroupSize of 5.

Results

Mitogenomic Variation

Sequencing reads from all ancient specimens showed the expected patterns of DNA fragmentation and deamination rates that were consistent with those of authentic aDNA (Supplementary Figure 1). We obtained 48 mitogenomes with at least 3-fold average coverage. We also obtained mitogenomes for 24 modern Orkney specimens (Supplementary Table 2). A total of 2135 SNPs (∼13% of mitogenome positions) were identified among all 545 samples – including the outgroup species Alaska pollock –: 1219 SNPs corresponded to informative sites and 916 SNPs were singletons (Supplementary Table 5).

Nucleotide diversity (π) between modern locations ranged between 0.002 and 0.003 (Table 1) and π of ancient samples did not vary from the values obtained in modern locations. The temporal comparison of specific sites (Quoygrew-Orkney and Haithabu-North East Arctic), showed limited significant differences between genetic statistics of temporally spaced ancient and modern locations (Supplementary Table 6 and Supplementary Figure 2), where Haithabu has significantly lower π and higher F compared to the North East Arctic (Supplementary Figure 2).

TABLE 1

Table 1. Estimates of genetic diversity statistics for Atlantic cod at whole mitogenomes from different marine locations or ecotypes in the North Atlantic (see Supplementary Tables 1, 2, 5).

Neutrality tests showed significant negative values for all Tajima’s D (TD) and F statistics in most locations, except for the western location Baffin Island, and the eastern locations Tvedestrand fjord and western Baltic (Table 1). Overall, there were 486 haplotypes —- including the outgroup—across all 545 samples, of which only 26 were shared between individuals (Figure 2 and Supplementary Table 7). Ancient CytB variation consisted of 7 different haplotypes, including four main haplotypes (A, C, D, and E) previously identified in modern mtDNA studies (Árnason, 2004; Jørgensen et al., 2018). Two novel variations of existing CytB haplotypes were found in western Baltic (haplotype ED) and North Sea (haplotype LI), while another 2 novel variations of existing CytB haplotypes were found among ancient samples (haplotypes LJ and TI). The most prevalent ancient haplotypes were A and E (∼ 40 and 38%, respectively, Supplementary Tables 1, 7), which were also commonly found in modern samples (Supplementary Table 8). The haplotype genealogy for whole mitogenome and CytB sequence data showed an extensive distribution of ancient samples across marine locations (Figure 2 and Supplementary Figure 4). Limited geographic mitogenome structure was observed, except for elevated divergence between western Atlantic and eastern Atlantic locations, and between locations in the western and eastern Baltic Sea and other eastern Atlantic locations (Figures 2B,C and Supplementary Figures 4B,C). A star-like topology is observed for the whole mitogenome and CytB genealogies (Figure 2 and Supplementary Figure 4).

FIGURE 2

Figure 2. Star-like haplotype genealogy based on complete mitogenomes of 545 samples of Atlantic cod. Only transversions (336) are depicted in this genealogy. (A) Haplotype genealogy colored by location. (B) Haplotype genealogy highlighting ancient specimens as a single group (in red; see Figure 1 and Supplementary Table 1; for specific distribution of sampling sites of ancient samples see Supplementary Figure 3). (C) Haplotype genealogy highlighting samples from the western (orange) and eastern (yellow) Baltic. (D) Haplotype genealogy highlighting western Atlantic (green) and Baffin Island (light green) specimens. Circle size is proportional to haplotype frequency. Black dots indicate intermediate haplotypes not found among the specimens and black lines connecting black circles correspond to mutation steps between haplotypes.

Demographic Patterns and Population Structure

The AMOVA assigned 7.58% of the variation between marine locations (including ancient samples as a single group) while 91.47% of the variation was represented between individuals (ΦCT = 0.076, p ≤ 0.001; ΦST = 0.085, p ≤ 0.000). Pairwise Φ_ST values (Figure 3 and Supplementary Tables 9, 10) showed significant differentiation levels between all ancient samples and western Atlantic, Baffin Island, western Baltic, eastern Baltic and Tvedestrand fjord. Ancient samples showed higher differentiation when compared to western Atlantic (Φ_ST = 0.117), and Baffin Island (Φ_ST = 0.192) in comparison to other eastern Atlantic locations. Among modern samples, western Atlantic, Baffin Island, western Baltic and eastern Baltic showed significant Φ_ST values when compared to all other locations (Supplementary Tables 9, 10). Φ_ST values were not significant between North Sea, North East Arctic, Norwegian coast and Ancient samples. Pairwise Φ_ST values between temporal spaced locations also showed no significant differentiation (Quoygrew and modern Orkney: Φ_ST = 0.000; p = 0.807; and Haithabu and North East Arctic: Φ_ST = 0.000; p = 0.456).

FIGURE 3

Figure 3. Heatmap of differentiation pairwise (Φ_ST) between 10 marine locations of Atlantic cod: western Atlantic (WA), Baffin Island (BI), western Baltic (WB), eastern Baltic (EB), Tvedestrand fjord (TV), North Sea (NS), North East Arctic (NE), Norwegian coast (NC), and Ancient specimens (AN). Dark brown colors indicate higher levels of differentiation (Φ_ST = 0.15–0.30). Light brown corresponds to lower levels of differentiation (Φ_ST = ≥ 0.10). Significant p-values (p ≤ 0.05) are indicated with a red *p-values are listed in Supplementary Tables 9, 10.

The time-calibrated Bayesian phylogeny for ancient and modern Atlantic cod samples resulted in 2 main clades with an estimated divergence from the most recent common ancestor at 220 kya (95% highest posterior density (HPD) = 194,780—249,980 kya; Figure 4). The first clade, which is not further divided, includes mitogenomes from 6 different widely scattered localities. The second clade was composed by 16 subclades with posterior probability > 0.8, with divergence times of ca. 100 kya. Clades and subclades in the phylogeny were not geographically structured, with the exception of most samples from western Atlantic, and most samples from western and eastern Baltic, which clustered together (Figures 2, 4).

FIGURE 4

Figure 4. Time calibrated collapsed Bayesian phylogeny of full mitogenomes from 525 Atlantic cod specimens using Alaska pollock (Gadus chalcogrammus) as an outgroup. Pie charts represent the marine locations distributed in each clade. Numbers beside pie charts indicate the number of individuals from each marine location distributed in each clade. Only branches with posterior probability > 0.8 are indicated next to the corresponding clade/subclade. For specific distribution of sampling sites of ancient samples see Supplementary Figure 6.

The Bayesian skyline analysis using different subsets of the data revealed a consistent pattern of step-wise population expansions followed by periods of constant population size. Expansions around 150, 50, and 10 kya are present in most subsets (Figure 5). A population expansion of Atlantic cod was identified ca. 50 kya in all subsets. The most recent expansion identified (around 10 kya), is only present in data sets that include clades with most Baltic Sea specimens (Figures 5A,B,E). Despite such differences, all analyses agree with a high and increasing female effective population size (Ne) of Atlantic cod (Ne = ca. 1,000,000—10,000,000) during the last ca. 100 kya, with highest estimates of Ne during the last few millennia (Figure 5).

FIGURE 5

Figure 5. Atlantic cod demographic history representing the fluctuations of effective female population size (Ne; in dark brown line) based on (A) data set I: all 524 sequences, (B) data set II: 476 modern sequences (excluding 48 ancient samples), (C) data set III: 273 sequences (excluding clades associated with most western Atlantic and Baltic Sea samples), (D) data set IV: 368 sequences (excluding the clade associated with most Baltic Sea samples), (E) data set V: 429 sequences (excluding clades associated with most western Atlantic samples) and (F) data set VI: 48 ancient sequences (excluding all modern samples). The 95% CI corresponds to the light brown area. The most recent 95% CI at time 0 is shown in each respective panel. Description of the clades excluded in data sets III, IV, and V can be found in Supplementary Figure 7.

Discussion

Here, we compared modern and ancient mtDNA diversity in Atlantic cod to investigate whether observed historical and contemporaneous census population declines (Hutchinson et al., 2003; Hylen et al., 2008; Limburg et al., 2008; Bartolino et al., 2012; Jonsson et al., 2016; Brattey et al., 2018) have had mitogenomic consequences. The temporal comparison of 48 ancient specimens to 496 modern conspecifics did not reveal consistent significant mitogenomic changes or measurable effective genetic population declines through time. Below, we discuss reasons why such genomic impacts may not be observed.

First, mitogenomic variation is high in modern Atlantic cod and is characterized by limited genetic differentiation between populations and incomplete lineage sorting over large spatial scales across its range in the North Atlantic (Jørgensen et al., 2018; Lait et al., 2018). Low observed genetic differentiation (Φ_ST) between Tvedestrand fjord and other Norwegian coastal locations, as well as between the North Sea, the North East Arctic and the Norwegian coast confirm this lack of geographic structuring over large parts of the eastern Atlantic (Figure 3). Indeed, the non-significant differentiation of all ancient samples with modern North Sea, North East Arctic and Norwegian coast is fully consistent with their presumed geographical origin and highlights the long-term lack of mtDNA structure in this region. Non-significant ΦST values between the Norwegian coastal locations and Tvedestrand fjord indicate possible recent migration of fish between such coastal communities and more restricted fjord populations (Knutsen et al., 2011). Compared to many terrestrial ecosystems, where populations can often be isolated by physical barriers—which restrain interbreeding and dispersal—(Hauser and Carvalho, 2008; Exadactylos et al., 2019), in marine ecosystems the absence of physical barriers promotes larger panmictic populations and Atlantic cod is no exception (Berg et al., 2016, 2017; Sodeland et al., 2016; Barth et al., 2017). Thus, a combination of low spatial resolution of mtDNA data as a result of continuous gene flow and connectivity may mask any local temporal erosion of mitogenomic diversity (Welch et al., 2012) in Atlantic cod.

Second, we determined high long-term estimates of effective population size (Ne = ca. 1,000,000–10,000,000; Figure 5), which is in agreement with earlier observations in Atlantic cod (Hardie et al., 2006; Therkildsen et al., 2010; Pinsky et al., 2021). Estimates of Ne can remain high in economically important fish species, even if their populations have experienced a large biomass decline (Hauser and Carvalho, 2008) since it takes hundreds of generations (i.e., depending on the generation time of the species; Amos and Balmford, 2001; Frankham et al., 2002) for the actual population numbers and breeding populations to be reflected in Ne (Hauser and Carvalho, 2008). In fact, simulations have shown that a population with theoretical Ne of 100 (which is several orders of magnitude lower than observed in Atlantic cod) would retain 75% of heterozygosity after 57 generations (Frankham et al., 2002; Welch et al., 2012). Given that such population declines take a very long time to lead to measurable genomic consequences, mtDNA—as a single locus—will have limited power to record such changes in populations of high Ne (Allentoft et al., 2014; Johnson et al., 2018; Thomas et al., 2019; Spencer, 2020). The absence of significant genetic changes in this study is consistent with the absence of such changes in genome-wide data using historical samples of Atlantic cod from the western and eastern Atlantic (Pinsky et al., 2021) and with the absence of such changes in mitogenomic data from other taxa that have similarly high estimates of Ne as Atlantic cod, such as the Pacific herring (Speller et al., 2012; Moss et al., 2016), the Hawaiian petrel (Welch et al., 2012) and even extinct species such as the New Zealand moa (Allentoft et al., 2014), the passenger pigeon (Murray et al., 2017) and the great auk (Thomas et al., 2019).

In contrast, temporal losses of mitogenomic diversity and/or declines in Ne have been reported in species that have suffered population fragmentation (e.g., resulting in small effective population sizes) or that have experienced limited connectivity, such as the steppe bison (Shapiro et al., 2004), the Scandinavian arctic fox (Nyström et al., 2006), cave bears (Stiller et al., 2010), the Iberian lynx (Casas-Marce et al., 2017), the Iberian salmon (Consuegra et al., 2002) and the common bream (Ciesielski and Makowiecki, 2005). Interestingly, a loss of haplotypic variation has been identified—using CytB sequence data—for a single period (i.e., 15th to 16th centuries, out of 6 temporal periods investigated) in an Icelandic population of Atlantic cod (Olafsdottir et al., 2014). There are two potential explanations for this discrepancy. First, nearly all substitutions that comprise the CytB haplotypes can be affected by post-mortem deamination (i.e., they consist of C > T and G > A substitutions). Most of the ancient sequences (90%) investigated in Olafsdottir et al. (2014) were obtained in a single round of PCR without evaluation of such post-mortem deamination. Therefore, such bias due to post-mortem damage cannot be excluded. Second, our sampling does not include many specimens from Iceland (Figure 1 and Supplementary Table 1), and it remains possible that—with 156 samples—a local effect has been observed in Olafsdottir et al. (2014), which we do not detect in our data.

Third, we do not observe major novel mtDNA lineages in the ancient data, nor observe a significant loss of such lineages over time. Instead, the majority of Atlantic cod mtDNA lineages observed in ancient and modern samples today have originated ca. 100–150 kya (Figure 4), during a period of population expansion (Figure 5). Therefore, the gain of such lineages—and associated population expansions—in Atlantic cod is more likely caused by changes in abundance driven by major historical climatic events such as eustatic oscillations in sea level, and the interglacial and warming periods experienced during the last glacial maximum ca. 23,000 kya (Bigg et al., 2008) and the Wisconsinan (ca. 110–120 kya) and Illinoian (ca. 200–130 kya) glaciations (Gibbard and Van Kolfschoten, 2005) as described by Lait et al. (2018). For instance, we only observe the most recent population expansion ca. 10 kya (Figures 5A,B,E) when including those mtDNA clades which are strongly associated with the Baltic Sea. The timing of this expansion is in agreement with the development of the Baltic Sea (ca. 7,000–8,000 years; Ojaveer et al., 2010; Wenne et al., 2020) which has led to genetically distinct Atlantic cod populations that have adapted to local environmental conditions (i.e., salinity and temperature; Johannesson and Andre, 2006; Berg et al., 2015; Wenne et al., 2020). Therefore, the observed changes in Ne reflect past population demography rather than recent and contemporary demographic changes (Lombal et al., 2020).

It is clear from zooarcheological evidence that Atlantic cod has periodically experienced intense exploitation in the distant past, particularly around the North Sea and the Baltic Sea (Barrett et al., 1999; Enghoff, 1999; Olson and Walther, 2007; Orton et al., 2011). This fishing pressure became even greater in the 19th and 20th centuries (e.g., Thurstan et al., 2010). Landings of Atlantic cod exceeded 4,000,000 tons during 1960–1990s in the North Atlantic Ocean (Shelton and Morgan, 2014). In particular, landings surpassed 600,000 tons in Iceland by ca. 1930s (Drinkwater, 2006), 354,000 tons in the North Sea during ca. 1970s (Bannister, 2004), 200–400,000 tons in the eastern Baltic during 1960–1990s (MacKenzie et al., 2002), 650,000 tons in North East Arctic between 1937 and 1938 up to 800–1,200,000 tons in ca. 1950s (Sætersdal and Hylen, 1964; Hylen, 2002). Such high levels of exploitation led to major reductions in present abundances of most Atlantic cod populations (i.e., Food and Agriculture Organization [Fao], 2020-2021a,b). Nonetheless, for the reasons discussed above, our results indicate that such population declines of Atlantic cod did not lead to a detectable impact on the mtDNA genome on the time scale we investigated here.

Taken together, our results highlight that historical and contemporaneous anthropogenic pressures such as commercial fisheries have had little impact on the ancient mitogenomic diversity of a wide-spread marine species with high gene flow such as Atlantic cod. Future ancient DNA studies should consider the inclusion of nuclear genomic data and extensive sampling on a local scale—considering a temporal comparison of specimens from the same geographical region—to assess the effects of climate and human exploitation with greater statistical power. Finally, our observations do not contradict evidence that overfishing has had negative consequences for the abundance of Atlantic cod and they do not oppose information about the important implications of genetic variation in evolutionary biology, ecology and conservation biology. Instead, our observations suggest that conservation management measures aimed toward the demographic recovery of Atlantic cod in the eastern Atlantic, if achievable by conservation management measures, will not be constrained by recent loss of historical mitogenomic variation.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ebi.ac.uk/ena, PRJEB42959.

Author Contributions

BS and JB: conceptualization, project design, and supervision. GF and LM-G: laboratory work and data curation. LM-G with input from GF, TO, LA, BS, and JB: formal analysis. HB and MS: analytical advice. JB, IY, IJ, JH, RN, DO, SH-D, US, AH, and RB: ancient Atlantic cod specimens. FB: modern Orkney specimens. JB, IY, IJ, JH, RN, DO, SH-D, US, AH, and RB: archeological context information. LM-G and BS: data visualization. JB, BS, KJ, SJ, and DO: funding acquisition. LM-G and BS with input from JB: writing—original draft. All authors writing—review and editing.

Funding

This work was supported by Research Council of Norway projects “Catching the Past” (262777) and “The Aqua Genome Project” (221734), Leverhulme Trust Project MRF-2013-065 and the European Union’s Horizon 2020 Research and Innovation Programme under the Marie Skłodowska-Curie grant agreement No. 813383.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Agata T. Gondek-Wyrozemska for help with processing the ancient specimens and Francis Neat for providing modern Orkney specimens. We also thank Oliver Kersten and Michael Matschiner for advice during analyses. Finally, we thank M. Skage, S. Kollias and A. Tooming-Klunderud at the Norwegian Sequencing Centre for sequencing and processing of samples. Analyses were performed on the SAGA Cluster using the resources and assistance from the SIGMA2 Metacenter, the Norwegian National Infrastructure for High Performance Computing and Data Storage.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2021.671281/full#supplementary-material

References

Allentoft, M. E., Heller, R., Oskam, C. L., Lorenzen, E. D., Hale, M. L., Gilbert, M. T. P., et al. (2014). Extinct New Zealand megafauna were not in decline before human colonization. Proc. Natl. Acad. Sci.U. S. A. 111, 4922–4927. doi: 10.1073/pnas.1314972111

PubMed Abstract | CrossRef Full Text | Google Scholar

Amos, W., and Balmford, A. (2001). When does conservation genetics matter? Heredity 87, 257–265. doi: 10.1046/j.1365-2540.2001.00940.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Árnason, E. (2004). Mitochondrial cytochrome b DNA variation in the high-fecundity Atlantic cod: trans-Atlantic clines and shallow gene genealogy. Genetics 166, 1871–1885. doi: 10.1093/genetics/166.4.1871