Historical Demographic Processes Dominate Genetic Variation in Ancient Atlantic Cod Mitogenomes

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.

Ancient DNA (aDNA) approaches have been successfully used to infer the longterm 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 mutationdrift 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 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 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(Barrett et al., , 2011Orton 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., 1998Sigurgí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  Table 1).
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.

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.

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 . 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 ) 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 alignmentincluding 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.
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 , 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 haplotypegenealogy graph using Fitchi (Matschiner, 2016) with the MLbased 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)  Finally, a Coalescent Bayesian Skyline (CBS) analysis was completed to reconstruct the demographic historyincluding 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.

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).
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).

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   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).
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).

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 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 lineagesand 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  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 1937and 1938up to 800-1,200,000 tons in ca. 1950s (Saetersdal and Hylen, 1964Hylen, 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.

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.