Ancient genomes of Sitka black-tailed deer show evidence for postglacial stepping-stone dispersal along the Pacific Northwest Coast of North America

doi:10.21203/rs.3.rs-5033480/v1

Download PDF

Research Article

Ancient genomes of Sitka black-tailed deer show evidence for postglacial stepping-stone dispersal along the Pacific Northwest Coast of North America

https://doi.org/10.21203/rs.3.rs-5033480/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Background:

The mule deer (Odocoileus hemionus) and its two distinct black-tailed deer (BTD) subspecies, Sitka and Columbian BTD, have a complex history in North America involving survival in Last Glacial Maximum (LGM) refugia, postglacial expansion along the Pacific Northwest Coast, evidence for incomplete lineage sorting and recent introgression between subspecies. Moreover, the differentiation process of the two black-tailed deer subspecies is poorly understood and could have been a consequence of the LGM. As such, they provide an exemplary system to explore patterns of population dynamics in response to climate change.

Results:

Here we analyzed genome-scale data from samples spanning the last 13,500 years to explore the evolutionary history of Sitka BTD in Southeast Alaska. Deer samples from Southeast Alaska older than 8,500 years ago shared a mitochondrial haplotype with mule deer, whereas samples younger than 6,000 years have the modern Sitka BTD haplotype. Discordantly, nuclear genomic data confirmed that all ancient individuals from Southeast Alaska are closely related to modern Sitka BTD, although the older group also shared ancestry with mule deer. Modern samples from Vancouver Island share more alleles with modern Sitka BTD than Columbian BTD. Our results support that they survived in the same glacial refugium south of the Cordilleran ice sheet, along today’s Oregon coast.

Conclusion:

The uneven deglaciation along the Northwest Pacific Coast following the LGM may have created temporary post-glacial refugia, or “stepping stones”, along the British Columbia Coast. Such dispersal, associated with genetic drift and isolation by distance, likely led to the emergence of the BTD subspecies, as well as the low genetic diversity observed in modern Sitka BTD.

paleogenomics

Southeast Alaska

Last Glacial Maximum

postglacial dispersal

During the Last Glacial Maximum (LGM), large parts of North America were covered by the Cordilleran and Laurentide Ice Sheets, reaching its maximum extent within the time frame 26 to 17 thousand years ago (ka) depending of the region [ka; 1, 2]. Ice-free areas, known as refugia, allowed species to survive during the glaciation. As the ice sheets retreated, migration routes became available, enabling species to recolonize previously glaciated regions and reconnect refugial populations. The initial proposed migration route that connected the north and south of the ice sheets involved a continental corridor between the Cordilleran and Laurentide Ice Sheets, spanning across Central Canada [3, 4]. However, recent studies suggest that while this corridor opened around 15 ka, it only became ecologically viable around 13.8 ka [3]. The northwestern Pacific coastal part of the Cordilleran Ice Sheet likely experienced an earlier and faster deglaciation, with areas such as the Alexander Archipelago in Southeast Alaska starting to become ice-free around 16 ka [2, 5, 6] and Vancouver Island around 18.5 ka [7]. This early deglaciation, combined with lower sea levels, could have created a viable coastal migration route thousands of years before the continental corridor.

Despite an earlier hypothesis suggesting ice-free conditions in areas of the Alexander Archipelago during a local LGM (lLGM; 20–17 ka), evidence from cosmogenic exposure dating has revealed that such areas that are above the modern sea level were indeed covered in ice during this period [2]. This finding is supported by genetic analyses conducted on ancient Southeast Alaska brown (Ursus arctos) and American black bear (U. americanus) mitogenomes, further corroborating complete ice coverage in these regions during the lLGM [8]. However, it is important to note that parts of the continental shelf, currently submerged, were exposed due to significantly lower sea levels. These submerged areas on the continental shelf might have remained ice-free during the local LGM.

The recovery of thousands of mammalian remains from multiple caves along the southern portion of the Alexander Archipelago, dating back to at least 45 ka, has provided valuable insights into the region's paleontological and climatic history [9]. The subfossil record includes evidence of humans and artifacts [10–12], as well as some of the oldest genetically confirmed New World dog remains [13]. This rich subfossil record reveals a clear transition in fauna, with certain species, such as brown and American black bears, persisting both before and after the local Last Glacial Maximum, while others, like caribou (Rangifer tarandus), became extirpated after the lLGM. Once deglaciation began, new species appeared, such as deer (Odocoileus hemionus). As a long-distance migrant herbivore [14], Odocoileus hemionus could have rapidly colonized these regions, making it an opportune species for understanding the timing of deglaciation along the northwestern Pacific Coast of North America. Additionally, it provides insights into the viability of the North Pacific Coastal route as ice retreated from the region and allows for the assessment of the impacts of climate change on large mammals in Southeast Alaska during the Late Pleistocene.

Odocoileus hemionus, a medium-sized cervid species, has a widespread distribution in Western North America, extending from Mexico to Alaska. The species is divided into five subspecies, further grouped into two main clusters: black-tailed deer (BTD) and mule deer (MD) [15]. Within the BTD group, there are two subspecies: O. h. columbianus (Columbian black-tailed deer; CBTD), spanning from northern California to the Central Coast of British Columbia, Canada, and O. h. sitkensis (Sitka black-tailed deer; SBTD) found along the Central Coast of British Columbia to Southeast Alaska, with recent (early 20th century) introductions to Haida Gwaii, British Columbia [16], and other coastal areas as far north as Kodiak Island, Alaska [17, 18]. The mule deer group encompasses three subspecies: O. h. hemionus (continental mule deer; MD), distributed across most of western North America, and O. h. cerrosensis and O. h. sheldoni that are endemic to Cedros Island and Tiburón Island, respectively, both in Mexico [15, 19].

Despite belonging to the same species, black-tailed and mule deer exhibit significant genetic divergence. Mitochondrial DNA loci reveal a genetic difference of 5–8% between the two main groups [20–23]. White-tailed deer (O. virginianus) and mule deer possess a similar mitochondrial haplotype which has been attributed to incomplete lineage sorting [24]. Nuclear analyses performed so far demonstrate less divergence among the taxa but notable differences persist [25–28].

According to the prevailing hypothesis, mule deer are believed to have survived the Last Glacial Maximum in multiple refugia located to the east of the Cascade Range, while black-tailed deer are thought to have survived in a single refugium, possibly on the western side of the Cascade Range, along today’s Oregon and Washington coast [20, 25]. The observed low genetic diversity in the BTD subspecies today has been hypothesized to be a consequence of the expansion out of a southern refugium until reaching its northern natural distribution, Southeast Alaska. The processes leading to today’s differentiation between the black-tailed deer subspecies are not fully understood. Hence, a northern coastal refugium, such as in the Alexander Archipelago in Southeast Alaska, has been suggested as a potential ice-free area during the Last Glacial Maximum for black-tailed deer [29] that could potentially explain the differentiation of the subspecies [15]. Black-tailed deer only appear in the subfossil record on the Vancouver Island and Haida Gwaii archipelago at approximately 13 calibrated thousands of years before the present (cal kyr B.P.) and in the Alexander Archipelago at around 9 cal kyr B.P. [9, 30, 31]. However, due to the incomplete subfossil record, the possibility of a northern refugium where deer went undetected cannot be ignored.

So far, all studies concerning Odocoileus hemionus evolution have been based on modern samples. Adding genetic data from ancient samples allows us to explore how glacial isolation in refugia may have affected the species, the colonization process following the ice sheet retreat, and changes in diversity over time, providing a deeper understanding of the species’ evolutionary history and its modern phylogeographic patterns. Here, we present results from genetic analyses of 21 ancient black-tailed deer from Southeast Alaska, Haida Gwaii, and Vancouver Island, spanning the last ~ 13.5 cal kyr B.P. We aim to explore the evolutionary history of Sitka black-tailed deer in the Alexander Archipelago, including whether the subspecies survived in the region during the Last Glacial Maximum, or if an ancestral black-tailed deer survived in southern refugia and expanded northwards. Understanding their out-of-refugium expansion can help provide insights into the opening and viability of a North Pacific coastal route and assess climate change impacts on herbivores during the Holocene.

We generated genomic data from 21 ancient black-tailed deer from Southeast Alaska, Haida Gwaii, and Vancouver Island, in addition to 13 new modern mule and black-tailed deer (Fig. 1, Table S1). All ancient samples were radiocarbon-dated, with calibrated ages ranging from ~ 0.8 cal kyr B.P. to ~ 13.5 cal kyr B.P. The ancient mitogenomes had an average sequencing depth ranging from 9.46x to 432x and were near-complete (> 90% coverage of width; Table S1). The ancient whole genomes had an average sequencing depth ranging from 0.01x to 1.65x and a breadth ranging from 1.6–71%, whereas modern samples had an average from 3.78x to 9.96x and a breadth ranging from 76–94% (Table S2). Ancient nuclear and mitochondrial genomes showed the expected pattern of degraded ancient DNA, with an increased rate of cytosine deamination at the 5’-end of the reads, and an increased rate of guanine to adenine substitutions close to the 3′- end of the reads, confirming the ancient DNA authenticity (Figure S1 and Figure S2).

Two distinct maternal lineages inhabited Southeast Alaska during the Holocene

Phylogenetic analyses of complete mitochondrial genomes confirmed previously reported groupings within modern Odocoileus, including two major clades, one encompassing black-tailed deer, and the other, mule deer plus white-tailed deer (Fig. 2a, see Figure S3 for the complete tree). One sample, MD5, possessed the mitochondrial haplotype of Columbian black-tailed deer, although whole genome data confirmed its taxonomic identity as mule deer. These two clades exhibited a deep split dating back to ~ 670 cal kyr B.P. (95% HPD ~ 478–930 cal kyr B.P.). Notably, ancient SE Alaska samples older than 8.5 cal kyr B.P. (referred to as old-SEAK) grouped with mule deer instead of black-tailed deer, while samples younger than 6 cal kyr B.P. (referred to as young-SEAK) grouped with modern Sitka black-tailed deer, suggesting a haplotype, and potentially subspecies, turnover between ~ 8.5 and 6 cal kyr B.P.

Within the black-tailed deer lineage, modern Columbian black-tailed deer from Oregon (CBTD2, CBTD3, CBTD4) formed a sister group to a highly supported monophyletic group comprising samples from Washington, British Columbia, and Alaska (including most ancient samples) that diverged ~ 28 cal kyr B.P. (95% HPD ~ 19.5–38 cal kyr B.P.) The most recent common ancestor of the Washington, British Columbia, and Alaska clade was estimated at ~ 22 cal kyr B.P. (95% HPD ~ 16.5–29 cal kyr B.P.). Despite low genetic diversity revealed by short branches, two main groups emerged inside the Washington, British Columbia, and Alaska clade: one consisting of ancient and modern samples from Vancouver Island, and the other comprising modern samples from Washington, as well as modern and ancient samples from British Columbia and Alaska. Modern Southeast Alaska samples were closely related to ancient samples younger than 6 cal kyr B.P. Two samples, CBTD_BC3 from the Central Coast of British Columbia, and SBTD4, from Haida Gwaii are placed as a sister group of a clade that only encompasses Sitka black-tailed deer from Alaska (modern and ancient younger than 6 cal kyr B.P.). Those two clades shared a last common ancestor with Sitka black-tailed deer from Alaska ~ 18.4 kyr B.P. (95% HPD ~ 12.6–26.2 cal kyr B.P.). Those two groups may suggest that Sitka black-tailed deer has two distinct maternal lineages, one from British Columbia and one from the Alexander Archipelago. All samples from Southeast Alaska (including those translocated to Kodiak Island) with a black-tailed deer haplotype shared a last common ancestor ~ 11.6 cal kyr B.P. (95% HPD ~ 8–16 cal kyr B.P.).

Maximum likelihood analyses showed a similar topology to the Bayesian tree within black-tailed deer, except for ancient samples from Haida Gwaii, which, in the Bayesian tree, are placed within a non-supported group sister to SE Alaska samples (Figure S4). In the maximum likelihood analysis, they were either sister to the clade encompassing samples from Washington, British Columbia, and SE Alaska, or sister to the Alaskan clade, receiving low bootstrap support in both cases.

The haplotype network revealed that the black-tailed deer and mule deer haplogroup differed by 624 substitutions, which correspond to a genetic difference of ~ 4.3% (Fig. 2b). Overall, the haplotype network captured similar relationships to the ones observed in the phylogenetic trees for O. hemionus, the main difference being the position of Columbian black-tailed deer (CBTD), and the samples from Vancouver Island (CBTD_BC). The southern group of Columbian black-tailed deer (CBTD) was now nested inside black-tailed deer. The SE Alaska group was separated from the ancient samples from Haida Gwaii cluster by at least 14 mutations and from the samples from the Central coastal British Columbia by at least 13 mutations.

The two distinct maternal lineages are genetically Sitka black-tailed deer.

To explore population structure from genome-wide deer data, genotype likelihoods from both modern and ancient samples were generated using ANGSD [32]. A dataset containing only O. hemionus samples was created and used for the principal component analysis [PCAngsd; 33]. Most modern samples clustered in three main groups, mule deer, Columbian black-tailed deer, and Sitka black-tailed deer (Fig. 3a). The majority of mule deer were placed into a tight cluster, suggesting limited genetic diversity within the subspecies. However, four samples (MD5, MD8, MD9, MD10) deviated towards the Columbian black-tailed deer cluster along both PC1 and PC2, hinting at potential historical and/or contemporary gene flow in areas where the subspecies come into contact [34].

Black-tailed deer displayed a non-perfect cline along PC1 and PC2, with modern Sitka black-tailed deer from SE Alaska at one extreme and Columbian black-tailed deer from Oregon at the other. In between these two clusters, there were samples from Haida Gwaii, Central British Columbia, Vancouver Island, and southern coastal British Columbia, in this order. Hence, the cline roughly follows the geographic distribution of black-tailed deer along the northwestern Pacific Coast, suggesting that the differentiation of the two black-tailed deer subspecies may a consequence of isolation by distance.

All ancient samples from Southeast Alaska clustered with modern Sitka black-tailed deer, confirming their taxonomic identity. However, ancient samples from the old-SEAK group (older than 8.5 cal kyr B.P.) exhibited a closer relationship to mule deer, supporting their mule deer mitochondrial haplotype. Ancient samples from Vancouver Island were positioned between the two black-tailed deer subspecies, relatively closer to modern samples from the same region. Those ancient samples have ages ranging from 13.2 cal kyr B.P. to 0.8 cal kyr B.P. and most of them are very closely related. To identify any possible trend related to age, a linear regression between PC1 and PC2 and the samples’ ages was performed. Regarding PC1, there was no correlation (r² = 0.14; Fig. 3b), while PC2 revealed a stronger trend (r² = 0.48; Fig. 3c), suggesting that older samples are more closely related to Sitka black-tailed deer than younger and modern samples from Vancouver Island.

Even though the variance explained by further principal components decreases (Figure S5), and no clear groups were observed, it is interesting to note that when analyzing PC2xPC3 and PC3xPC4, most modern black-tailed deer were still separated by subspecies. However, ancient samples seemed to cluster by groups and genome coverages on PC2xPC3. Specifically, SEAK3 and SEAK4 grouped together and separated from the other SEAK samples, while VI2, VI5, and VI9 placed into a different cluster, separated from the other ancient samples from Vancouver Island. Those samples exhibited higher sequencing coverage (ranging from 0.54x to 1.5x) when compared with samples from the same regions (ranging from 0.01x to 0.2x). PC3-PC6 linear regressions showed a similar trend to PC2, with younger samples closer to Columbian black-tailed deer; however, the r² was low (0.11 to 0.23; Figure S6)

Signatures of mitochondrial-nuclear discordance

Single nucleotide polymorphisms (SNPs) from modern individuals were called using ANGSD [32], excluding ancient individuals due to their low coverage (depth ranging from 0.01 to 1.65x). A distance-based phylogenetic network was constructed using the NeighborNet method [Figure 4a; 35]. The two black-tailed deer and mule deer groups formed an interconnected network, characterized by extra edges among branches. This complexity in the network suggests data conflict, potentially arising from incomplete lineage sorting (ILS), admixture, or recent divergence. Sitka black-tailed deer (SBTD) and black-tailed deer from Vancouver Island and Central British Columbia (CBTD_BC) branched out of the black-tailed deer network, with Sitka black-tailed deer nested inside this group. The maximum-likelihood tree recovered a similar relationship as the network (Figure S7).

To further investigate and visualize incongruent phylogenomic signals in the genomic data we next employed Twisst [Topology Weighting by Iterative Sampling of Subtrees; 36]. We used non-overlapping windows of 50 SNPs, resulting in a total of 5296 windows. We separated modern samples into five groups that were identified based on the Principal Component Analyses (Fig. 3a). White-tailed deer (WTD) was used as an outgroup. In the three top topologies, the CBTD_BC samples grouped with SBTD, even though they were harvested within the Columbian black-tailed deer geographic distribution area (Fig. 4b). Notably, in three out of five topologies with the highest weights, CBTD_BC grouped with SBTD, only in the fifth-highest weighted tree did CBTD_BC and CBTD group together. This suggests that CBTD_BC individuals are closer related to Sitka black-tailed deer than Columbian black-tailed deer.

Because incongruent phylogenetic signals can be caused by introgression or incomplete lineage sorting, to identify the source of the incongruence, we next performed QuiBL [Quantifying Introgression via Branch Lengths; 37]. We generated 286 phylogenies based on 1,000 SNPs in non-overlapping windows. We found that the proportion of discordant trees that arose from introgression was 14.28%, which corresponds to 3.33% of all phylogenies, indicating that most of the discordance is caused by incomplete lineage sorting. Introgression was observed mostly inside the white-tailed, black-tailed, and mule deer lineages. Some gene flow was detected between black-tailed and mule deer; however, no gene flow was detected between the two Odocoileus species (Figure S8).

Results from exploring signatures of gene flow in further detail

Ancestry sharing – To estimate the ancestry of the studied samples, NGSAdmix [38], a model-based clustering analysis based on genotype likelihoods, was employed. The optimal number of ancestral groups (K) was determined as K = 4 using CLUMPAK [39]. Two distinct ancestral groups were identified for mule deer (MD), one with a higher representation in the southern United States and Mexico, and a second more concentrated in Canada (Fig. 5a, see Figure S9 for K = 3 to K = 7). Columbian black-tailed deer from Oregon and Washington (CBTD) had a distinct ancestral group, different from Sitka black-tailed deer (SBTD). The older samples from SE Alaska shared ancestry with the northern mule deer group (~ 15%) and modern Sitka black-tailed deer. The younger samples belonged to the same ancestral group as modern Sitka black-tailed deer (Fig. 5b). Interestingly, modern samples from Vancouver Island (CBTD_BC1 and 2) shared ancestry with both black-tailed deer subspecies in an approximately equal amount (Fig. 5a), while ancient samples from Vancouver Island had a higher proportion of the Sitka black-tailed deer ancestral group (Fig. 5b). The correlation between Sitka ancestry percentage and age of ancient Vancouver Island samples returned a r² = 0.42 (Fig. 5c), suggesting a decrease in Sitka ancestry proportion over time, consistent with the observed trend in PCA.

f3 – To test for definitive evidence of admixture, we performed the f3 statistic (Fig. 6a-d; Figure S10). The difference between ancient and modern samples from Vancouver Island and Southeast Alaska was notable. While modern samples from Vancouver Island (CBTD_BC1 and CBTD_BC2; Fig. 6a) did not exhibit any significant value (Z score below − 3) once they were tested as a target, the three ancient samples with a depth of coverage above 0.5x (VI2, VI5, and V9; Fig. 6b) were highly significant whenever a black-tailed deer was one of the sources. However, if a CBTD individual is one of the sources, significant values were only observed if the second source was a deer from Vancouver Island (modern and ancient), Central Coast of British Columbia, and Alaska (modern and ancient). It is noteworthy that highly significant results were obtained when ancient samples were placed as targete, with mule and black-tailed deer (except for CBTD individuals) used as sources. Even though, modern samples from Vancouver Island do not return significant values, a similar trend can be observed with less positive values.

Modern samples from Southeast Alaska just yield a significant Z score if another Sitka black-tailed deer is one of the sources (Fig. 6c). However, ancient samples from Southeast Alaska (Fig. 6d) exhibited a similar pattern to ancient samples from Vancouver Island than with modern Sitka black-tailed deer. The main difference between the ancient groups was when a Columbian black-tailed deer was placed as one of the sources, as it returned no significant values. This may be due to the isolation since the Last Glacial Maximum.

When mule deer was a source and the sample was from around where black-tailed and mule deer geographic ranges meet (MD5, MD8, MD9, MD10), the Z values were significantly negative if any black-tailed deer was a target, indicating gene flow among them. All other mule deer had low values when two mule deer were sources, but highly positive values if a black-tailed deer was a source instead.

D-statistics – We next evaluated the amount of drift shared among modern and ancient black-tailed deer populations along the Northwestern Pacific Coast using D-statistics, by computing D (CBTD, SBTD; X, WTD). Where X were samples from Vancouver Island (modern and ancient), Southeast Alaska (ancient). All tested combinations returned significantly positive D and Z scores, indicating gene flow with Sitka black-tailed deer (Fig. 6e). Hence, samples from Vancouver Island and Southeast Alaska shared more alleles with modern Sitka black-tailed deer than with Columbian black-tailed deer.

Although the levels of genetic drift shared between modern Sitka black-tailed deer and the different ancient samples from Vancouver Island were similar, they exhibited some variation. The linear regression of the D score and the age of those samples returned an r² value of 0.69, indicating a strong correlation, with older samples sharing a higher amount of drift with Sitka black-tailed deer than the younger ones (Fig. 6f). The old-SEAK group from Southeast Alaska (SEAK3, SEAK4, and SEAK5) displayed similar values among each other and with ancient Vancouver Island samples, suggesting that those two regions were colonized from a similar refugia. On the other hand, there was an observable increase in the amount of drift shared between modern Sitka black-tailed deer and the young-SEAK (SEAK6, SEAK7, and SEAK8) compared to the previously mentioned groups.

TreeMix – A maximum likelihood drift tree was generated to infer multiple splits and mixture patterns among modern Odocoileus hemionus individuals, using TreeMix [40]. The tree without any migration edge, was similar to the maximum likelihood phylogeny, explaining most of the variance in the dataset (Figure DS4). However, 0.26% of the variance was not captured by this tree, and based on the addition of 1 to 10 migration events and OptM results [41], the optimal number of migrations was M = 1. This first admixture edge, which was found consistently throughout all TreeMix runs, indicates an admixture event from black-tailed deer (CBTD_BC3) to a mule deer (MD8), samples that exhibited ancestry sharing with the three O. hemionus subspecies (Fig. 5a). Other admixture events included edges inside black-tailed deer and mule deer lineages, as well as from black-tailed deer to mule deer (Figure S11).

∆-statistics – To further explore gene flow within Odocoileus, we performed the recently proposed ∆-statistics [42] to test for the direction of introgression among modern samples, testing the asymmetric tree, A = ((((SBTD, CBTD_BC), CBTD), MD), WTD) (Fig. 7a). Approximately 60% of the 100 kb windows support the species tree (Fig. 7b), while in analyses based on larger windows (500 kb), only 36% of windows show no sign of introgression (Fig. 7c). When analyzing scenarios of gene flow between different subspecies, the direction from black-tailed deer to mule deer is predominant. The scenario CBTD_BC ◊ mule deer (2 ◊ 4) was supported by ~ 8% of the 100 kb windows, while the opposite scenario was supported by only ~ 2%. A change in the proportion of these two scenarios is noticeable when looking at 500 kb windows, which may favor more recent events: 17% of windows support introgression from black-tailed deer to mule deer, while the opposite direction is supported by ~ 9%. Gene flow from Sitka black-tailed deer to mule deer (1◊ 4) shows similar percentages of windows with both window sizes (~ 2.4%). However, there is a slight difference as the opposite scenario received slightly more support with 500 kb windows. For both window sizes, significant bidirectional scenarios are also observed within the black-tailed deer group. A very small number of windows support introgression between black-tailed and white-tailed deer. These admixture events corroborate the results from TreeMix analyses (Figure S11).

Levels of genetic diversity

To look for changes in genetic diversity over time, we calculated the heterozygosity of modern and ancient black-tailed deer, as well as mule deer. Ancient samples from the Alexander Archipelago had a similar heterozygosity; however, the old group (older than 8.5 cal kyr B.P.) had a wider breadth when compared to the younger group (Fig. 8). It is possible to observe a sharp decline in heterozygosity when comparing ancient and modern samples from Southeast Alaska. Vancouver Island samples older than 10 cal kyr B.P. had a slightly higher heterozygosity than samples younger than 5.5 cal kyr B.P. Both of those groups had values higher than modern individuals from the same region. Hence, there was a loss in genetic diversity over time in the Alexander Archipelago and on Vancouver Island, and it could have been associated with repeated founder events, drift, and their isolation on islands for a long period. Mule deer and Columbian black-tailed (CBTD) deer had comparable values, and higher than their insular conspecifics. Sitka and Columbian black-tailed deer exhibited less divergence (F_ST = 0.11) than Sitka black-tailed and mule deer (F_ST = 0.23).

F_ST outliers

Finally, we investigated whether genomic regions of differentiation may help identify functional divergence between Sitka and Columbian black-tailed deer. Genomic regions were considered outliers when their fixation index (Fst) was above the 99% percentile (0.6). No clear islands of divergence between the two subspecies were found (Figure S12). Gene ontology analysis results indicated that most of the enriched regions are involved in biological processes and genes, including sensory perception (e.g. OR7C2, OR7A17), reproduction (e.g. DMC1), immune response (e.g. LCK), development process (e.g. DOCK7, SNORC, TMEM223, TR2B), and cancer-related genes (e.g. SOX7, CDH19) (Table S4).

Although the two black-tailed deer subspecies are genetically distinct, their evolutionary history is not yet fully understood. For example, previous studies of modern Odocoileus hemionus based on mitochondrial DNA and microsatellite data supported a hypothesis that both black-tailed deer subspecies survived the Last Glacial Maximum along today’s Oregon and Washington coast, and dispersed north once conditions became favorable [20, 25]. However, those previous studies were solely based on modern individuals, which only represent a portion of past diversity. As such they may incompletely describe the evolutionary history of species as multiple demographic and evolutionary events can result in the same genetic signature in modern samples. For example, the low genetic diversity observed in modern Sitka black-tailed deer could have been a consequence of a bottleneck due to refugia survival in small populations, expansion out of refugia, or both [43]. Hence, analyzing ancient samples may significantly improve the understanding of the evolution of black-tailed deer. The earliest presence of black-tailed deer in the fossil record in Vancouver Island [~ 13.5 cal kyr B.P.; 30], Haida Gwaii [~ 12.8 cal kyr B.P.; 31], and the Alexander Archipelago in Southeast Alaska [~ 9.2 cal kyr B.P.; 9] corroborate a hypothesis of LGM occupation in southern refugia, followed by a northward dispersal along the Northwest Pacific coast. However, an incomplete fossil record may have left deer undetected from those regions. Due to the considerably lower sea level following the LGM in comparison with today, refugia along outer coastal areas of British Columbia or Southeast Alaska may have existed in areas that are now submerged [2, 44], and survival in separate LGM refugia could potentially explain the distinctiveness of the two black-tailed deer subspecies.

Deer in British Columbia

Black-tailed deer have inhabited Vancouver Island since at least 13.5 cal kyr B.P. However, there is a notable gap in the fossil record from 13.5 to 5.5 cal kyr B.P., with only one sample recovered so far during this time period (VI10; ~10.4 cal kyr B.P.). Despite this gap, the close genetic relatedness of ancient samples across this time gap suggests that black-tailed deer may have been present in the region during this period, yet are unsampled. In this case, deer may have survived the Younger Dryas, a cooling period from ~ 12.8 to 11.5 cal kyr B.P. [45], in the region. Interestingly, ancient samples from Vancouver Island appear to have had higher genetic diversity than modern individuals, indicating a loss of genetic diversity over time, probably due to long-term island isolation and potentially population contractions during unfavorable periods in the Holocene.

Previous research has identified that extant Columbian black-tailed deer from Vancouver Island and Gabriola Islands belong to a different ancestral group when compared to mainland deer [25]. The island group is closely related to ancient individuals from Vancouver Island and shares a significant amount of alleles with modern Sitka black-tailed deer suggesting a possible origin through northward dispersal. During population expansions, individuals at the leading edge of dispersal are more successful in passing their alleles to future generations [46]. This phenomenon, known as gene surfing, increases the likelihood of certain alleles becoming fixed in the population. Ancestral black-tailed deer alleles may have become fixed in deer populations inhabiting the region, and due to geographic isolation, these alleles persist in the deer populations of Vancouver Island today. This process could explain the genetic distinctiveness observed in the island group compared to mainland Columbian black-tailed deer.

Our results from principal component analysis, D-statistics, and ADMIXTURE analyses show a clear trend in linear regression, where older samples from Vancouver Island are closer related to Sitka black-tailed deer than younger and modern samples are (modern VI samples share about 50% of their ancestry with each black-tailed deer subspecies). The increased sharing of VI alleles over time with Columbian black-tailed deer (CBTD) may be due to secondary dispersal waves.

All ancient samples from Haida Gwaii possess a black-tailed deer mitochondrial haplotype; however, their placement in the Bayesian and maximum likelihood phylogenetic trees is incongruent. In the Bayesian analyses, all three ancient samples from Haida Gwaii are placed in a monophyletic clade sister to a clade that encompasses modern samples from Haida Gwaii, Central coastal British Columbia, and ancient and modern samples from Alaska, whereas in the maximum likelihood tree, they were either sister to a clade that encompasses modern samples from Haida Gwaii, Central coastal British Columbia, and ancient and modern samples from Alaska, or sister to the Vancouver Island clade. These incongruences may indicate that Haida Gwaii samples were only isolated for a short period of time and did not differentiate from other black-tailed deer, or incomplete fossil record and gaps in sampling of modern inviduals.

The fossil record indicates that deer had a very short presence in Haida Gwaii. Around 13.5 cal kyr B.P. due to lower sea level, Haida Gwaii was a much larger island stretching eastwards towards the mainland. This facilitated deer to occupy the archipelago [Figure 9c; 31, 47]. Deer disappeared from the Haida Gwaii fossil record ~ 12.8 cal kyr B.P., which coincides with the beginning of the Younger Dryas period [~ 12.8 to 11.5 cal kyr B.P.; 45]. The extirpation of deer in Haida Gwaii also coincides with the extirpation of brown bears [31]. Based on the fossil record, deer never recolonized the Haida Gwaii, until the early 20th century, when a population was introduced from the British Columbia mainland coast, and it is still present there today [16].

The arrival of deer in Southeast Alaska

Southeast Alaska marks Sitka black-tailed deer’s northernmost native distribution (individuals from the Alexander Archipelago were introduced to the Kodiak Archipelago in the early 20th century [17, 29]), but the timing of its arrival is still an open question. All deer samples from Southeast Alaska younger than 6 cal kyr B.P. shared a last common ancestral Sitka black-tailed deer matriline around 12 cal kyr B.P. which could suggest the time interval that deer arrived in Southeast Alaska and became isolated in the Alexander Archipelago. However, because older samples from SE Alaska (~ 9.2–8.5 cal kyr B.P.) possessed a mule deer mitochondrial haplotype, this divergence timing is likely underestimated.

Based on paleoclimate records and climate models, Praetorius, et al. [48] suggested periods when conditions could have been favorable for human dispersal from Beringia to the south of the ice sheets along the Northwestern Pacific Coast. These periods, likely also coinciding with favorable conditions for deer dispersals from south of the ice sheet to Southeast Alaska, occurred from 24.5 to 22 cal kyr B.P., 16.4 to 14.8 cal kyr B.P., and 13 to 11.7 cal kyr B.P. Considering the estimated divergence time of maternal haplotypes, the entrance of deer in the fossil records from Vancouver Island, Haida Gwaii, and Southeast Alaska, and the history of sea levels in the Alexander Archipelago between 13.5 (Fig. 9c) and 11.5 cal kyr B.P. (Fig. 9d), it is probable that deer could have reached the Alexander Archipelago during the last period, 13 to 11.7 cal kyr B.P. Following this period, the sea level rose up to 20 meters higher than today on the outer coast [44], potentially hindering dispersal to the Archipelago.

Two distinct Sitka black-tailed deer matrilineal lineages have inhabited SE Alaska. Ancient samples dated from ~ 9.2–8.5 cal kyr B.P. share their mitochondrial genome with mule deer, whereas all samples younger than 6 cal kyr B.P. possess a modern Sitka black-tailed deer mitochondrial haplotype. Although all studied samples before 8.5 cal kyr B.P. have a mule deer haplotype, due to an incomplete fossil record and a small sample size, deer with the Sitka black-tailed deer haplotype may have also been present in the Alexander Archipelago at this time. Nevertheless, Heaton and Grady [9] reported that one of the oldest O. hemionus recovered in Southeast Alaska (SEAK3), had an unusual antler with unique protuberances at its base, and it was different from modern Sitka black-tailed deer antlers observed in the region. Moreover, the antlers did not bifurcate, they had a sinuous above the eye sockets wider than any modern Sitka black-tailed deer, and the lower jawbone was thinner than in modern individuals. It is possible that the older SE Alaska deer group had a different morphology and was part of an isolated population, as suggested by Heaton and Grady [9].

Today, Sitka black-tailed deer in the Alexander Archipelago are the Odocoileus hemionus subspecies with the lowest genetic diversity. Modern Sitka black-tailed deer possess a significantly lower heterozygosity when compared to ancient individuals, although the “old” and “young” ancient groups yield similar heterozygosity values. The individual from Haida Gwaii in our study (SBTD4) had a higher heterozygosity than extant samples from Alexander Archipelago and Kodiak Island, indicating different levels of genetic diversity within modern Sitka black-tailed deer, possibly between island and mainland populations. The lower values observed in the Alexander Archipelago and Kodiak Archipelago (introduced from the Alexander Archipelago) may be due to long-term island isolation, whereas deer in Haida Gwaii were introduced from the mainland stock and have been isolated for a shorter period, allowing them to have a higher genetic heterozigozity when compared individulas from Alaska.

Introgression has been reported among black-tailed deer and mule deer, but has particularly focused on Columbian black-tailed and mule deer [34, 49, 50]. By contrast, little discussion has centered on modern introgression between Sitka black-tailed and mule deer. This may stem from limited geographic contact between these subspecies due to the Coastal Mountains acting as a significant barrier. However, our results with ∆-statistics were able to capture introgression from Sitka black-tailed (and the CBTD_BC group) to mule deer. Even though samples from Vancouver Island did not possess mule deer alleles or the mule deer haplotype, potentially due to an incomplete fossil record and small sample size, the mule deer mitochondrial haplotype may have been present in Vancouver Island individuals at the time of glacial retreat and could have become fixed by drift once they became isolated in the Alexander Archipelago. In this case, the introgression with mule deer could have been an older event that affected ancestral black-tailed deer individuals and diminished over time. However, our QuiBL results showed that most of the incongruence present in our modern dataset arose due to incomplete lineage sorting. Hence, it is likely that the presence of mule deer ancestry in modern populations of mule deer and black-tailed deer, as well as in ancient samples from Southeast Alaska is a consequence of a combination of past introgression between those lineages and incomplete lineage sorting.

There is a gap of deer in the fossil record from 8.5 to 6 cal kyr B.P. in Southeast Alaska, after which a matrilineal lineage turnover may have occurred, and a different lineage of Sitka black-tailed deer appeared. It is important to note that fossils from other species have been recovered from this gap, such as black bears and otters [9], suggesting that an incomplete fossil record may not be the only explanation for this gap. Furthermore, around the same period, brown bears disappeared from the fossil record on the southern islands [8, 9]. Wilcox, et al. [51] reconstructed the paleoclimate in the region during the last 13.5 cal kyr B.P. using δ¹⁸O concentrations in speleothems from caves in the Alexander Archipelago. δ¹⁸O is used as a proxy to understand past climate changes, as higher concentrations indicate colder periods whereas lower concentrations indicate warmer periods. The transition from the Bølling-Allerød warming period, around 12.9 to 12.7 ka, to the Younger Dryas period was marked by a rapid increase in δ¹⁸O concentration, which indicates a decrease in temperature. The δ¹⁸O remained constant until around ~ 9.5 cal kyr B.P. when a relatively rapid drop in δ¹⁸O was reported, suggesting an increase in temperature. However, at 8.5 ka, the δ¹⁸O rapidly increased, indicating a rapid cooling event. It is important to note that this was the most abrupt change in δ¹⁸O since 12.5 cal kyr B.P. [51]. This rapid cooling could have affected the deer population inhabiting the area, causing a bottleneck or a complete extirpation of deer from the region, followed by later recolonization. Southeast Alaska marks the northernmost native distribution of Sitka black-tailed deer, and even today, winter remains an important limiting factor for the species in the region [52, 53].

A postglacial northwards stepping-stone dispersal of black-tailed deer

Although prevous studies identified genes potentially involved with the functional divergence of white-tailed and mule deer [27] and mule and black-tailed-deer [54] it was suggested that selection may have played a relatively small role in deer speciation [27]. We were not able to identify candidate genes that may have been associated with an adaptive divergence of the Sitka and Columbian black-tailed deer lineages. Some enrichment of genes related to sensory perception and immune response was found, but such genes are generally found enriched in mammals and vertebrates [55, 56]. Signatures of positive selection may be confounded by modern gene flow and the recent split of Columbian and Sitka black-tailed deer that most likely occurred soon after the Last Glacial Maximum. Therefore, future scans of genes associated with functional differentiation between the two subspecies warrant more detailed analysis including a denser sampling of populations, including introgressed individuals between the two subspecies. For example, we just compared Columbian black-tailed deer from Oregon and denser sampling along the Northwest Pacific Coast is needed to fully understand the differentiation of the two black-tailed deer subspecies.

Modern Columbian black-tailed deer from Oregon, Washington, and southern British Columbia (CBTD) possess the highest genetic diversity of black-tailed deer, followed by Vancouver Island individuals, and Southeast Alaska. The decrease in genetic diversity along the Northwest Pacific Coast supports the hypothesis that an ancestral population of the two subspecies survived in a single refugium, likely around today’s Oregon and Washington coasts. This is coupled with the fact that deer from Vancouver Island, Central Coast of British Columbia (CBTD_BC), and Alaska (SBTD) are closely related to each other, and the timing of deer entrance in the fossil record corroborates this hypothesis. Once deglaciation began, and conditions became favorable, deer may have dispersed northwards in one main migration wave. Based on a mitochondrial molecular clock, this initial dispersal could have happened as early as ~ 21 cal kyr B.P., which is comparable to dates obtained by previous studies [19 cal kyr B.P.; 20]. However, the presence of the Cordilleran Ice Sheet along the coast during this time may have impacted their dispersal northwards.

The Cordilleran Ice Sheet reached its maximum southwestern extent (Puget Lowland, Washington) around 17 cal kyr B.P., which was followed by a fast retreat [57, 58]. However, coastal areas along outer Vancouver Island and parts of the Haida Gwaii Archipelago were likely already ice-free by 18.5–17 cal kyr B.P. [Figure 9a,b; 7, 31, 47], whereas deglaciation only began in Southeast Alaska around 16 − 15 cal kyr B.P. [2, 6]. Such temporally asymmetric deglaciation along the Northwest Pacific coast might have delayed black-tailed deer in their initial postglacial migration wave moving north as the ice sheet could have acted as a barrier during the dispersal. Until approximately 17 cal kyr B.P. (Fig. 9b), the Cordilleran Ice Sheet is thought to have covered most of the Northwest Pacific Coast; however, the lower sea level at the beginning of deglaciation on the outer coast beginning at least as early as 14,500 cal kyr B.P. may have provided ecologically viable areas along the coast for deer and other mammalian species in this period. This unique deglaciation pattern may have created temporary ice-free refuges surrounded by unsuitable habitat and ice, similar to a “stepping-stone” landscape [59]. Hence, the initial dispersal northward along the Northwest Pacific Coast may be characterized as stepping-stone migration, with temporary genetic isolation along the British Columbia coast until black-tailed deer reached Southeast Alaska. Due to gene surfing during population expansion, ancestral alleles present in black-tailed deer that were at the front edge of the dispersal wave could have been fixed by genetic drift during isolation. Such a stepping-stone style of dispersal, associated with expansion out of a postglacial refugium, may explain the cline observed with PCA and the decrease of genetic diversity as a result of northward dispersal along the coast, as well as the higher allele sharing between deer from British Columbia (Vancouver Island and Central British Columbia) and modern Sitka black-tailed deer.

The Sitka black-tailed deer subspecies lost significant genetic diversity compared to its close modern and past relatives, challenging the reconstruction of their evolutionary history solely based on modern samples. Ancient samples from Southeast Alaska and British Columbia have provided important insights, not only clarifying the complex evolutionary history of Odocoileus hemionus in North America but also shedding further light on the opening and viability of the coastal route along the Northwest Pacific Coast.

Our findings do not support the survival of deer in Southeast Alaska during the LGM. Instead, our data provide evidence that a black-tailed deer ancestor survived in a single refugium along today’s Oregon coast. As deer dispersed northward, ice sheets persisted in some areas along the British Columbia coast, delaying their expansion and creating temporary post-glacial refuges along the Canadian coast. Such stepping-stone dispersal associated with genetic isolation and drift in these refuges contributed to reduced genetic diversity in modern deer inhabiting Southeast Alaska.

Sampling and radiocarbon dating

Twenty-one subfossil deer samples from the Alexander Archipelago, Haida Gwaii, and Vancouver Island were analyzed (Table S1; Fig. 1). All subfossils had been previously identified based on morphology, with most of the samples from Southeast Alaska classified as deer without a subspecies assignment, and two samples identified as deer or caribou [9]. Of the three samples from Haida Gwaii, two were morphologically identified as ungulate, and one as deer [31]. All Vancouver Island samples had been morphologically identified as deer [30]. We also generated whole genome data from six modern black-tailed deer from Alaska, British Columbia, and Oregon, and seven modern mule deer from different parts of its North American distribution (Table S2; Fig. 1). In addition to the newly generated data, we downloaded complete mitogenome data and whole genome data available in public databases (Table S2).

All ancient samples analyzed in this study had been previously ¹⁴C radiocarbon-dated based on bone collagen [9, 30, 31]. The radiocarbon dates were calibrated using the IntCal20 calibration curve in OXCAL v4.3 [60], reporting 2-sigma (Table S1).

DNA extraction, PCR amplification, mitochondrial enrichment, and shotgun sequencing

DNA from subfossils was extracted in a dedicated cleanroom facility at the University at Buffalo, separated from any processing of modern samples. Ancient DNA extractions followed the protocol described in Dabney, et al. [61] and the modifications in da Silva Coelho, et al. [13]. DNA from modern samples was extracted using either a salt extraction procedure [62, 63] or a DNeasy Blood and Tissue kit (Qiagen).

To obtain an initial genetic-based taxonomic identification of the samples, PCR reactions were performed, either adding 21µL H₂O, 0.4 µM of each forward and reverse primer, and 2 µL of extracted genomic DNA to each GE Illustra PuReTaq Ready-To-Go PCR bead (GE Healthcare) or by adding 2.5 µl of 10 × PCR buffer (Applied Biosystems, USA), 2.5 mM of each dNTP (Applied Biosystems), 25 mM of MgCl₂ (Applied Biosystems), 5–10 U µl^− 1 of AmpliTaq Gold DNA polymerase (Applied Biosystems), 10 µM each of the forward and reverse primers, 2 µl of the genomic DNA and 17.4 µl of H₂O. Because all the samples had been identified as deer/ungulate, we designed taxon-specific primers and amplified two regions of mitochondrial DNA cytochrome b (using primers UB285F 5’-ATCACAAATCTCCTCTCAGC-3’/UB286R 5’-TGGACTATAGCRAGTGCTGCGA-3’ and UB287F 5’-CCAGCAAACCCACTCAAYAC-3’/UB288R 5’-ACTCCTTAGTTTATTTGGA-3’, respectively), and one region of the control region (UB289F 5’-TACCATCCCTGAAACCA-3’/UB290R 5’-TGGCCCTGAAGAAAGAACCAG-3’, respectively). PCR products were purified using a MinElute PCR Purification kit (Qiagen), and Sanger sequenced directly using the same primers as in the PCR reaction.

To assemble complete mitochondrial genomes of ancient samples, library preparation and mitochondrial DNA target enrichment were performed by Arbor Daicel Biosciences (http://www.arborbiosci.com). DNA from twelve ancient samples was target-enriched using a bait panel developed from white-tailed deer (NCBI accession NC-015247). Four of the sample libraries followed the standard MYBAITS v. 3.0 protocol, with equal masses of libraries pooled, bead-templated, and sequenced on the Ion Proton platform. Following sequencing, reads were demultiplexed, quality-trimmed, and filtered using the default settings on the Ion Torrent suite v. 4.4.3. Following sequencing, reads were demultiplexed, quality-trimmed and filtered using the default settings on the Ion Torrent suite v. 4.4.3. For the remaining eight samples, Truseq dual-barcoded libraries were prepared without sonication, using the blunt-end ligation module from the NEBNext Fast DNA library preparation kit (New England BioLabs) with an extended double-time treatment and blunt-end adapters synthesized by Arbor Biosciences and paired-end sequenced on an Illumina HiSeqX platform.

Additionally, low-depth Illumina shotgun sequencing was performed for 15 ancient deer samples (six samples from Southeast Alaska and nine from Vancouver Island), in addition to 13 modern O. hemionus individuals. Single-stranded libraries from ancient samples were prepared using ssDNA2 [64], and double-stranded libraries from modern samples were prepared using the NEBNext Fast DNA library preparation kit (New England BioLabs). All libraries were sequenced on an Illumina HiSeqX platform. One sample (SBTD5) was sequenced on an Illumina MiSeq following the manufacturer's recommended protocols, conducted at the U.S. Geological Survey's Alaska Science Center, Anchorage, Alaska.

Genome mapping assembly and DNA degradation assessment

Due to the substantial mitochondrial divergence between mule and black-tailed deer, coupled with the absence of publicly available Sitka black-tailed deer mitogenomes when initiating this study, we assembled a de novo mitochondrial genome from a Sitka black-tailed deer harvested from Southeast Alaska (SBTD5). This assembly was achieved using Novoplast [65], a tool designed for de novo assembly of organellar genomes. To assemble mitogenomes from ancient samples, Illumina adapters were initially trimmed with AdapterRemoval v. 2.3.2 [66], reads shorter than 20 bp were discarded, and two base pairs were trimmed from the 5’ and 3’ ends. Ancient samples were aligned separately against reference mitogenomes (our own black-tailed deer mitogenome and mule deer NC-020729.1). Reads were aligned to the reference using the bwa-aln algorithm v. 0.7.17 [67], with maximum edited distance (-n) set to 0.01, the maximum number of gap open (-o) to 2, and the seed (-l) to 16500, unmapped reads were extracted using samtools v. 1.16.1 [68] and mapped with BWA-mem v. 0.7.17 [69] using default settings. PCR duplicates were removed using the MarkDuplicates tool in Picard v. 2.25.0 (http://broadinstitute.github.io/picard/) using the lenient validation stringency. The consensus was called using samtools mpileup [70] and the default settings. Mapping statistics were calculated with BEDTools v. 2.30.0 [71]. Reads from modern samples were mapped together using BWA-mem [69] under the default settings. PCR duplicate removal, consensus calling, and mapping statistics were performed following the same pipeline described above. To assemble the nuclear genome, the same pipeline as above was used, with the reads aligned against the white-tailed deer reference genome Ovir.te_1.0 [GCF_002102435.1; 72].

To analyze the damage pattern and assess DNA authenticity of ancient samples, we repeated the same mapping pipeline for ancient samples as described above, however, we omitted the trimming step of two base pairs from the 5’ and 3’ ends. We used mapDamage2 v. 2.0.8 [73], which uses an approximate Bayesian estimation of the damage patterns, to assess DNA degradation patterns of the reference-mapped assemblies (mitochondrial and nuclear genomes) of each ancient deer individual.

Mitochondrial genome analyses

In addition to the 34 newly assembled mitogenomes for this study, we downloaded complete mitogenomes of modern deer species from the tribe Odocoileini: Mazama sp. (18), Pudu sp. (2), Blastocerus dichotomous (1), Ozotoceros bezoarctos (2), Hippocamelus antisensis (1), Odocoileus virginianus (9), O. hemionus hemionus (2) from the National Center for Biotechnology Information (NCBI) Genbank database (Figure S3). Reads from an additional two O. h. sitkensis, four O. h. columbianus, ten O. h. hemionus, and 18 O. virginianus were downloaded from the NCBI Sequence Reads Archive (SRA; Table S1), and the mitogenomes were assembled following the same pipeline described above for modern samples. The total 106 sequences were aligned in MAFFT [74], with a manual inspection in GENEIOUS v. 2023.0.4 to remove tandem repeats from the control region. To exclude unaligned regions, we performed G-block in SeaView v. 5.0.5 [75]. Maximum likelihood analyses were performed using IQ-TREE v. 2.2.1 [76] with 1000 bootstraps using the GTR + G + I model that was chosen as the best model by IQ-TREE. To obtain estimated divergence data, BEAST v. 2.7.1 [77] was performed, using the GTR + G + I substitution model and a constant-size coalescent model, trees were sampled every 1000 states from a total of 70 million states, age calibration of divergence time estimation was performed by adding tip dates using calibrated radiocarbon dates from ancient samples. South American deer (Mazama sp.) were used as an outgroup. An additional dataset composed only of Odocoileus hemionus was used to generate a TCS parsimony network using POPART [78].

Whole genome analyses

In addition to the genome data from 15 ancient and 13 modern samples generated in this study, 32 genomes, including from O. virginianus, were downloaded from public repositories (Table S2). Three distinct datasets were created for downstream genomic analyses (Table S3). The first included all modern and three ancient samples with an average sequencing depth above 0.5x (DS1; N = 48), the second included ancient and modern O. hemionus only (DS2; N = 43), and the third dataset only included modern samples (DS3; N = 45). Quality filters, genotype likelihood, and SNPs call were performed using ANGSD v. 0.94 [32]. For all datasets, the following quality filters were employed: adjustment for excessive mismatches (-C 50), probability of a base pair being misaligned (-baq 1), reads with multiple best hits were removed (-uniqueOnly = 1) and reads with a flag above 255 were removed (remove_bads = 1). To calculate genotype likelihood and call SNPs, the following parameters were used: estimation of the posterior genotype probability based on the allele frequency as prior (doPost = 1), genotype likelihood obtained with GATK algorithm (GL = 2), p-value for a site considered a SNP (SNP_pval = 1e-6), the frequency at which of major and minor allele were estimated (-domaf 1), major and minor alleles were outputted (-doGeno = 1), only scaffolds over 1 Mb were used (-rf), SNPs were called and output into a BCF file (-doBCF). For DS2, the minimum depth of the genotypes was set to 2x, and max 500x (-setMinDepthInd and -setMaxDepthInd). For DS2, each site needed to be present in at least 32 individuals. The BCF file was converted to VCF using BCFtools v. 1.14 [79]. For DS3, only SNPs with a minimum depth of 3x, and max depth of 500x were kept using VCFtools v. 0.1.16 [minDP, and maxDP; 80]. Missing data were removed in VCFtools (--max-missing 1.0).

Analyses of population structure and phylogenetic reconstruction

To visualize population structure among modern and ancient O. hemionus, we performed principal component analyses (PCA) using PCAngsd v. 0.982 [33] under the default settings with the dataset DS2.

A maximum likelihood tree was constructed with DS3 using IQ-TREE v. 2.2.1 [76], 1000 bootstraps, and ascertainment bias correction, and the TVM + F + ASC + G4 model was chosen as the best model by IQ-tree. To visualize incongruences in the dataset, a distance-based phylogenetic network using the NeighborNet method was created using SplitsTree4 v. 4.19.2 [35, 81]. Incongruences are shown as extra edges that can indicate incomplete linage sorting and/or admixture [82]. Distances were corrected using LogDet.

To investigate incongruences in more detail, we employed Twisst [Topology Weighting by Iterative Sampling of Subtrees; 36], a method that constructs alternative phylogenies based on sliding windows of SNPs across the genome and quantifies their contribution to the complete tree. We utilized non-overlapping windows of 50 SNPs. Maximum likelihood trees for each window were constructed using RAxML v. 8.0 [83] and the raxml_sliding_windows.py script (https://github.com/simonhmartin/genomics_general). Modern samples were categorized into five groups based on Principal Component Analyses: Columbian black-tailed deer from Washington, Oregon, and southern British Columbia (CBTD), samples from Vancouver Island and Central British Columbia (CBTD_BC), mule deer (MD), and Sitka black-tailed deer (SBTD), with white-tailed deer (WTD) used as an outgroup. Incongruences in phylogenies can be caused by incomplete lineage sorting and introgressions, and to determine the source of these signals, we also employed QuiBL [Quantifying Introgression via Branch Lengths; 37]. QuiBL analyzes trees generated from non-overlapping sliding windows and estimates the proportion of introgressed loci by analyzing independent triplets based on branch lengths. We utilized non-overlapping windows of 1 kb and generated maximum likelihood trees as described above.

Tests of admixture among Odocoileus hemionus groups

NGSAdmix

To estimate ancestry sharing within ancient and modern O. hemionus (dataset DS2), we conducted a cluster analysis using NGSadmix [38] with K ranging from 3 to 7, where K represents the number of ancestral sources. To determine the optimal number of ancestral populations, we executed 10 replicates for each K and analyzed the results using the Evanno [84] method through Clumpak [39].

f3-statistics

To test for more definitive evidence of admixture, we next performed the f3 statistic using AdmixTools v. 7.0.1 [85]. We computed f3(C; A, B), where C is the individual being tested (target), and A and B are the sources. A Z value below − 3 indicates that the target is admixed or has a close genetic relationship (population sharing) with one of the sources due to the presence of segments that are identical by descent [86]. We calculated f3 values from all combinations of modern samples and ancient samples with a depth above 0.5x based on SNPs generated with ANGSD (both transversions and transitions were included). Z values were converted into p-values and corrected using p.ajdust and FDR methods [87, 88] and then reconverted to Z values.

TreeMix

A maximum likelihood drift tree was constructed using TreeMix v. 1.13 [40], with only modern individuals, using white-tailed deer as an outgroup. TreeMix evaluates historical mixtures between populations based on genome-wide allele frequency data. TreeMix was run with migration edges ranging from 0 to 10 and the –noss flag activated, to disable correction for small sample size, and prevent overcorrection. To identify the optimal number of migration events, 10 replicates were executed for each migration value. The output from these replicates was then utilized as input for OptM [41], which employs the Evanno method [84] to assess and determine the most suitable number of migration events.

∆-statistics

To explore the direction of introgression events during the evolutionary history of Odocoileus, we employed ∆-statistics [42]. A VCF-file only including modern individuals was recoded to a variant-major additive component file (--recode A-transpose) using PLINK version 1.9 [89]. Using site-frequency data, we tested the asymmetric combination ((((SBTD, CBTD_BC), CBTD), MD), WTD) using non-overlapping windows (correction off) of 100 kb and 500 kb.

D-statistics (ABBA-BABA)

To further estimate the extent of gene flow, we performed a four-sample test of admixture, D-statistic (ABBA-BABA), based on genotype likelihoods using ANGSD v. 0.94 [-doAbbababa2; 32, 90]. All modern and ancient individuals were included, transitions and transversions were used, the block size was set to the default (5 Mb), and only including scaffolds over 1 Mb were. Unadmixed white-tailed deer was used as an outgroup, based on Kessler, et al. [27]. To test whether the ancient samples and modern samples from British Columbia (CBTD_BC) are more closely related to Sitka or Columbian black-tailed deer we computed D (CBTD, SBTD; X, WTD) for individuals, where CBTD represents a Columbian black-tailed deer individual from Oregon, Washington and southern British Columbia, SBTD a modern Sitka black-tailed deer, WTD a white-tailed deer individual, and X the sample that is being tested. A positive result indicates gene flow between the sample tested and Sitka black-tailed deer, while a negative result indicates gene flow with Columbian black-tailed deer. As there was no significant difference in D, regardless of which CBTD or SBTD sample was used they were combined into two groups, one per subspecies.

Genetic diversity

To estimate the genetic diversity of modern and ancient Odocoileus hemionus, we used –het in VCFtools [80]. Individual heterozygosity was calculated as “(N_SITE - O(HOM)) /N_SITE” from the output. We also calculated the pairwise fixation index (F_ST) between both modern black-tailed deer subspecies and mule deer to measure population differentiation using the VCFtools [80] options "--weir-fst-pop" and "--fst-window-size 10000".

F_ST outliers

To examine putative islands of divergence of Sitka black-tailed deer from its conspecific, Columbian black-tailed deer, we performed a genome-wide F_ST outlier analysis. Windows with an F_ST above the 99th percentile were treated as outliers. To identify genes in the outlier windows, we compared the window coordinates with the genome annotation of the white-tailed deer provided by DNAZoo [72], using BEDtools intersect v. 2.3 [71]. Genes were extracted and blasted against the cow genome (UAR2.0; GCA_002263795.4) to create a study file. To look for significantly enriched regions across their genome, we used GOATOOLS [91].

Deglaciation and exposure of continental shelves along the Northwest Pacific Coast

To determine periods during the Late Pleistocene that could have been viable for deer migration based on ice sheet retreat, bathymetric data, and sea level fluctuation, we downloaded the North American Deglaciation Isochrones (NADI-1) database from 19, 17, 13.5, and 11.5 cal kyr B.P. from Dalton, et al. [47], and incorporated bathymetric data from the Northwest Pacific Coast [92], and sea level history from the Alexander Archipelago [44].

Conflict of interest statement

The authors declare no conflicts of interest.

Supplementary Tables

Table S1.

Voucher information and sequencing statistics for deer subfossil samples analyzed in this study.

Table S2.

Sample information for modern deer assembled for this study.

Table S3

Nuclear SNP data sets for analyses.

Table S4

Gene ontology analysis results showing all significantly enriched regions (p_bonfferoni < 0.05)

Author Contribution

C.L designed the study; F.A.S.C., C.M.T., D.M., D.F., J.H., E.L., S.T., J.B, and T.H.H. provided samples and/or generated the genomic data; K.K.P., D.M., and J.B the geological and glacial context of the ice-sheet along the Northwest Pacific Coast; T.H.H, J.B., and D.M. provided the cave and/or paleontological context. F.A.S.C and C.L. analyzed the data; F.A.S.C and C.L. wrote the manuscript with contributions from all authors.

Acknowledgment

The authors thank the University of Alaska Museum Earth Sciences collection (UAMES) and the Museum of Southwestern Biology (MSB) for loan of the specimens. The authors are also grateful to Tongass National Forest archaeologists Jane Smith, Gina Esposito, and Jackie de Montigny, University of South Dakota student field assistants Frank Andy Klock, Nathan Carter, Brandon Silver, Louis Rezac, Clarissa Ford, Alex Santos, and Christy Heaton. This research was supported with funding from the National Science Foundation (DEB award #1556565, EAR award #1854550, 9870343, and 0208247, and EF award #2221988). Unpublished genome assemblies and sequencing data for Ovir.te_1.0 (GCF_002102435.1) was used with permission from the DNA Zoo Consortium (dnazoo.org). The Hakai Institute supported work conducted on Vancouver Island that provided the ancient deer samples reported on here. Gwaii Haanas National Park Reserve supported work conducted on Haida Gwaii that provided ancient deer samples from Haida Gwaii.

Data Availability

The mitochondrial genome sequences generated in this study are deposited in the NCBI GenBank database with Accession no. XXXXXXX to XXXXXXX. Raw reads for ancient and modern deer are deposited in the NCBI Sequence Read Archive with Accession no. XXXXXX.

Clark PU, Dyke AS, Shakun JD, Carlson AE, Clark J, Wohlfarth B, Mitrovica JX, Hostetler SW, McCabe AM. The last glacial maximum. Science. 2009;325:710–4.
Walcott CK, Briner JP, Baichtal JF, Lesnek AJ, Licciardi JM. Cosmogenic ages indicate no MIS 2 refugia in the Alexander Archipelago, Alaska. Geochronology. 2022;4:191–211.
Clark J, Carlson AE, Reyes AV, Carlson EC, Guillaume L, Milne GA, Tarasov L, Caffee M, Wilcken K, Rood DH. The age of the opening of the Ice-Free Corridor and implications for the peopling of the Americas. Proceedings of the National Academy of Sciences 2022, 119:e2118558119.
Heintzman PD, Froese D, Ives JW, Soares AER, Zazula GD, Letts B, Andrews TD, Driver JC, Hall E, Hare PG, Jass CN, MacKay G, Southon JR, Stiller M, Woywitka R, Suchard MA, Shapiro B. Bison phylogeography constrains dispersal and viability of the Ice Free Corridor in western Canada. Proceedings of the National Academy of Sciences 2016, 113:8057–8063.
Lesnek, Briner JP, Lindqvist C, Baichtal JF, Heaton TH. Deglaciation of the Pacific coastal corridor directly preceded the human colonization of the Americas. Sci Adv. 2018;4:eaar5040.
Lesnek AJ, Briner JP, Baichtal JF, Lyles AS. New constraints on the last deglaciation of the Cordilleran Ice Sheet in coastal Southeast Alaska. Quatern Res. 2020;96:140–60.
Hebda CFG, McLaren D, Mackie Q, Fedje D, Pedersen MW, Willerslev E, Brown KJ, Hebda RJ. Late Pleistocene palaeoenvironments and a possible glacial refugium on northern Vancouver Island, Canada: Evidence for the viability of early human settlement on the northwest coast of North America. Q Sci Rev. 2022;279:107388.
da Silva Coelho FA, Gill S, Tomlin CM, Papavassiliou M, Farley SD, Cook JA, Sonsthagen SA, Sage GK, Heaton TH, Talbot SL. Ancient bears provide insights into Pleistocene ice age refugia in Southeast Alaska. Molecular Ecology; 2023.
Heaton TH, Grady F. The late Wisconsin vertebrate history of Prince of Wales Island, southeast Alaska. Ice age cave faunas North Am. 2003;2:17.
Aqil A, Gill S, Gokcumen O, Malhi RS, Reese EA, Smith JL, Heaton TT, Lindqvist C. A paleogenome from a Holocene individual supports genetic continuity in Southeast Alaska. Iscience 2023, 26.
Dixon EJ. Late Pleistocene colonization of North America from Northeast Asia: New insights from large-scale paleogeographic reconstructions. Mobility and ancient society in Asia and the Americas. Springer; 2015. pp. 169–84.
Lindo J, Achilli A, Perego UA, Archer D, Valdiosera C, Petzelt B, Mitchell J, Worl R, Dixon EJ, Fifield TE. Ancient individuals from the North American Northwest Coast reveal 10,000 years of regional genetic continuity. Proceedings of the National Academy of Sciences 2017, 114:4093–4098.
da Silva Coelho FA, Gill S, Tomlin CM, Heaton TH, Lindqvist C. An early dog from southeast Alaska supports a coastal route for the first dog migration into the Americas. Proceedings of the Royal Society B 2021, 288:20203103.
Robinette WL. Mule deer home range and dispersal in Utah. J Wildl Manag 1966:335–49.
Heffelfinger JR, Latch EK. Origin, classification, and distribution. Ecology and Management of Black-tailed and Mule Deer of North America. CRC; 2023. pp. 3–24.
Burgess BT, Irvine RL, Russello MA. Population genomics of Sitka black-tailed deer supports invasive species management and ecological restoration on islands. Commun Biology. 2022;5:223.
Smith RB. History and current status of Sitka black-tailed deer in the Kodiak Archipelago. In Sitka black-tailed deer, proceedings of a conference. 1979: 184–195.
Paul TW. Game transplants in Alaska. Alaska Department of Fish and Game, Division of Wildlife Conservation; 2009.
Latch EK, Heffelfinger JR. Genetics informs meaningful intraspecific taxonomy: the black-tailed and mule deer complex. Anim Prod Sci. 2022;63:1615–22.
Latch EK, Heffelfinger JR, Fike JA, Rhodes OE Jr. Species-wide phylogeography of North American mule deer (Odocoileus hemionus): cryptic glacial refugia and postglacial recolonization. Mol Ecol. 2009;18:1730–45.
Cronin MA, Vyse ER, Cameron DG. Genetic relationships between mule deer and white-tailed deer in Montana. J Wildl Manag 1988:320–8.
Cronin MA. Mitochondrial-DNA phylogeny of deer (Cervidae). J Mammal. 1991;72:553–66.
Carr SM, Hughes GA. Direction of introgressive hybridization between species of North American deer (Odocoileus) as inferred from mitochondrial-cytochrome-b sequences. J Mammal. 1993;74:331–42.
Klicka LB, Najar N, Vázquez-Miranda H, Zink RM. Relationships among North American deer based on mitochondrial DNA and ultraconserved elements, with comments on mito-nuclear discordance. Mammal Res 2024.
Latch EK, Reding DM, Heffelfinger JR, Alcalá-Galván CH, Rhodes OE Jr. Range‐wide analysis of genetic structure in a widespread, highly mobile species (Odocoileus hemionus) reveals the importance of historical biogeography. Mol Ecol. 2014;23:3171–90.
Cronin MA. Mitochondrial and nuclear genetic relationships of deer (Odocoileus spp.) in western North America. Can J Zool. 1991;69:1270–9.
Kessler C, Wootton E, Shafer ABA. Speciation without gene-flow in hybridizing deer. Mol Ecol. 2023;32:1117–32.
Kessler C, Shafer ABA. Genomic Analyses Capture the Human-Induced Demographic Collapse and Recovery in a Wide-Ranging Cervid. Mol Biol Evol 2024, 41.
Latch E, Amann R, Jacobson J, Rhodes O Jr. Competing hypotheses for the etiology of cryptorchidism in Sitka black-tailed deer: an evaluation of evolutionary alternatives. Anim Conserv. 2008;11:234–46.
McLaren D, Wigen R, Fedje D, Dyck A, Hebda CF, Morien E, Pedersen MW, Willerslev E, Rutledge LY, Barrera MA. Late Pleistocene Faunal Assemblages from Karst Cave Settings on Northern Vancouver Island, Canada. PaleoAmerica. 2023;9:216–36.
Fedje D, Mackie Q, McLaren D, Wigen B, Southon J. Karst caves in Haida Gwaii: Archaeology and paleontology at the Pleistocene-Holocene transition. Q Sci Rev. 2021;272:107221.
Korneliussen TS, Albrechtsen A, Nielsen R. ANGSD: analysis of next generation sequencing data. BMC Bioinformatics. 2014;15:1–13.
Meisner J, Albrechtsen A. Inferring population structure and admixture proportions in low-depth NGS data. Genetics. 2018;210:719–31.
Latch EK, Kierepka EM, Heffelfinger JR, RHODES OE JR. Hybrid swarm between divergent lineages of mule deer (Odocoileus hemionus). Mol Ecol. 2011;20:5265–79.
Bryant D, Moulton V. Neighbor-net: an agglomerative method for the construction of phylogenetic networks. Mol Biol Evol. 2004;21:255–65.
Martin SH, Van Belleghem SM. Exploring Evolutionary Relationships Across the Genome Using Topology Weighting. Genetics. 2017;206:429–38.
Edelman NB, Frandsen PB, Miyagi M, Clavijo B, Davey J, Dikow RB, García-Accinelli G, Van Belleghem SM, Patterson N, Neafsey DE, Challis R, Kumar S, Moreira GRP, Salazar C, Chouteau M, Counterman BA, Papa R, Blaxter M, Reed RD, Dasmahapatra KK, Kronforst M, Joron M, Jiggins CD, McMillan WO, Di Palma F, Blumberg AJ, Wakeley J, Jaffe D, Mallet J. Genomic architecture and introgression shape a butterfly radiation. Science. 2019;366:594–9.
Skotte L, Korneliussen TS, Albrechtsen A. Estimating individual admixture proportions from next generation sequencing data. Genetics. 2013;195:693–702.
Kopelman NM, Mayzel J, Jakobsson M, Rosenberg NA, Mayrose I. Clumpak: a program for identifying clustering modes and packaging population structure inferences across K. Mol Ecol Resour. 2015;15:1179–91.
Pickrell J, Pritchard J. Inference of population splits and mixtures from genome-wide allele frequency data. Nat Precedings 2012:1–1.
Fitak RR. OptM: estimating the optimal number of migration edges on population trees using Treemix. Biology Methods Protocols. 2021;6:bpab017.
Leppälä K, Coelho FAS, Richter M, Albert VA, Lindqvist C. Five-leaf generalizations of the < em > D-statistic reveal the directionality of admixture. bioRxiv 2024:2024.2002.2024.581856..
Loog L, Thalmann O, Sinding MHS, Schuenemann VJ, Perri A, Germonpré M, Bocherens H, Witt KE, Samaniego Castruita JA, Velasco MS. Ancient DNA suggests modern wolves trace their origin to a Late Pleistocene expansion from Beringia. Mol Ecol. 2020;29:1596–610.
Baichtal JF, Lesnek AJ, Carlson RJ, Schmuck NS, Smith JL, Landwehr DJ, Briner JP. Late Pleistocene and early Holocene sea-level history and glacial retreat interpreted from shell-bearing marine deposits of southeastern Alaska, USA. Geosphere. 2021;17:1590–615.
Kaufman DS, Axford YL, Henderson AC, McKay NP, Oswald WW, Saenger C, Anderson RS, Bailey HL, Clegg B, Gajewski K. Holocene climate changes in eastern Beringia (NW North America)–A systematic review of multi-proxy evidence. Q Sci Rev. 2016;147:312–39.
Waters JM, Fraser CI, Hewitt GM. Founder takes all: density-dependent processes structure biodiversity. Trends Ecol Evol. 2013;28:78–85.
Dalton AS, Dulfer HE, Margold M, Heyman J, Clague JJ, Froese DG, Gauthier MS, Hughes AL, Jennings CE, Norris SL. Deglaciation of the north American ice sheet complex in calendar years based on a comprehensive database of chronological data: NADI-1. Q Sci Rev. 2023;321:108345.
Praetorius SK, Alder JR, Condron A, Mix AC, Walczak MH, Caissie BE, Erlandson JM. Ice and ocean constraints on early human migrations into North America along the Pacific coast. Proceedings of the National Academy of Sciences 2023, 120:e2208738120.
Jackson HH. A hybrid deer of the F 2 generation. J Mammal 1921:140–3.
Haines ML, Luikart G, Amish SJ, Smith S, Latch EK. Evidence for adaptive introgression of exons across a hybrid swarm in deer. BMC Evol Biol. 2019;19:1–17.
Wilcox PS, Spötl C, Honkonen J, Edwards RL. A Walker switch mechanism driving millennial-scale climate variability. Innov Geoscience. 2023;1:100026.
Hanley TA. Relationships between Sitka black-tailed deer and their habitat. US Department of Agriculture, Forest Service, Pacific Northwest Forest and &#8230.
Jones G, Bunnell F. Response of black-tailed deer to winters of different severity on northern Vancouver Island. In Proceedings of a symposium on fish and wildlife relationships in old-growth forests, 12–15 April 1982. Juneau, American Institute of Fishery Research Biologist Morchead City … 5–396.
Powell JH, Amish SJ, Haynes GD, Luikart G, Latch EK. Candidate adaptive genes associated with lineage divergence: identifying SNPs via next-generation targeted resequencing in mule deer (Odocoileus hemionus). Mol Ecol Resour. 2016;16:1165–72.
Nielsen R, Bustamante C, Clark AG, Glanowski S, Sackton TB, Hubisz MJ, Fledel-Alon A, Tanenbaum DM, Civello D, White TJ, Sninsky J, Adams J, Cargill MD. A Scan for Positively Selected Genes in the Genomes of Humans and Chimpanzees. PLoS Biol. 2005;3:e170.
Bernatchez L, Landry C. MHC studies in nonmodel vertebrates: what have we learned about natural selection in 15 years? J Evol Biol. 2003;16:363–77.
Steffen ML. New age constraints for human entry into the Americas on the north Pacific coast. Sci Rep. 2024;14:4291.
Porter SC, Swanson TW. Radiocarbon age constraints on rates of advance and retreat of the Puget lobe of the Cordilleran ice sheet during the last glaciation. Quatern Res. 1998;50:205–13.
Simberloff D, Farr JA, Cox J, Mehlman DW. Movement Corridors: Conservation Bargains or Poor Investments? Conserv Biol. 1992;6:493–504.
Ramsey CB, Lee S. Recent and planned developments of the program OxCal. Radiocarbon. 2013;55:720–30.
Dabney J, Knapp M, Glocke I, Gansauge M-T, Weihmann A, Nickel B, Valdiosera C, García N, Pääbo S, Arsuaga J-L, Meyer M. Complete mitochondrial genome sequence of a Middle Pleistocene cave bear reconstructed from ultrashort DNA fragments. Proceedings of the National Academy of Sciences 2013, 110:15758–15763.
Medrano JF, Aasen E, Sharrow L. DNA extraction from nucleated red blood cells. Biotechniques 1990, 8.
Sonsthagen SA, Talbot SL, White CM. Gene flow and genetic characterization of Northern Goshawks breeding in Utah. Condor. 2004;106:826–36.
Gansauge M-T, Gerber T, Glocke I, Korlević P, Lippik L, Nagel S, Riehl LM, Schmidt A, Meyer M. Single-stranded DNA library preparation from highly degraded DNA using T4 DNA ligase. Nucleic Acids Res. 2017;45:e79–79.
Dierckxsens N, Mardulyn P, Smits G. NOVOPlasty: de novo assembly of organelle genomes from whole genome data. Nucleic Acids Res. 2017;45:e18.
Schubert M, Lindgreen S, Orlando L. AdapterRemoval v2: rapid adapter trimming, identification, and read merging. BMC Res Notes. 2016;9:1–7.
Li H, Durbin R. Fast and accurate short read alignment with Burrows–Wheeler transform. bioinformatics 2009, 25:1754–1760.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078–9.
Li H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv preprint arXiv:13033997 2013.
Li H, Durbin R. Inference of human population history from whole genome sequence of a single individual. Nature. 2011;475:493.
Quinlan AR. BEDTools: the Swiss-army tool for genome feature analysis. Current Protocols in Bioinformatics 2014, 47:11.12. 11-11.12. 34.
Dudchenko O, Batra SS, Omer AD, Nyquist SK, Hoeger M, Durand NC, Shamim MS, Machol I, Lander ES, Aiden AP, Aiden EL. De novo assembly of the < i > Aedes aegypti genome using Hi-C yields chromosome-length scaffolds. Science. 2017;356:92–5.
Jónsson H, Ginolhac A, Schubert M, Johnson PL, Orlando L. mapDamage2.0: fast approximate Bayesian estimates of ancient DNA damage parameters. Bioinformatics. 2013;29:1682–4.
Katoh K, Rozewicki J, Yamada KD. MAFFT online service: multiple sequence alignment, interactive sequence choice and visualization. Brief Bioinform. 2019;20:1160–6.
Gouy M, Guindon S, Gascuel O. SeaView version 4: a multiplatform graphical user interface for sequence alignment and phylogenetic tree building. Mol Biol Evol. 2010;27:221–4.
Nguyen L-T, Schmidt HA, Von Haeseler A, Minh BQ. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol Biol Evol. 2015;32:268–74.
Bouckaert R, Vaughan TG, Barido-Sottani J, Duchêne S, Fourment M, Gavryushkina A, Heled J, Jones G, Kühnert D. De Maio N: BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis. PLoS Comput Biol. 2019;15:e1006650.
Leigh JW, Bryant D. POPART: full-feature software for haplotype network construction. Methods Ecol Evol. 2015;6:1110–6.
Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, Whitwham A, Keane T, McCarthy SA, Davies RM. Twelve years of SAMtools and BCFtools. Gigascience. 2021;10:giab008.
Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, Handsaker RE, Lunter G, Marth GT, Sherry ST. The variant call format and VCFtools. Bioinformatics. 2011;27:2156–8.
Huson DH, Bryant D. Application of phylogenetic networks in evolutionary studies. Mol Biol Evol. 2006;23:254–67.
Burgon JD, Vences M, Steinfartz S, Bogaerts S, Bonato L, Donaire-Barroso D, Martínez-Solano I, Velo-Antón G, Vieites DR, Mable BK. Phylogenomic inference of species and subspecies diversity in the Palearctic salamander genus Salamandra. Mol Phylogenet Evol. 2021;157:107063.
Stamatakis A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics. 2014;30:1312–3.
Evanno G, Regnaut S, Goudet J. Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol. 2005;14:2611–20.
Patterson N, Moorjani P, Luo Y, Mallick S, Rohland N, Zhan Y, Genschoreck T, Webster T, Reich D. Ancient admixture in human history. Genetics. 2012;192:1065–93.
Lan T, Leppälä K, Tomlin C, Talbot SL, Sage GK, Farley SD, Shideler RT, Bachmann L, Wiig Ø, Albert VA. Insights into bear evolution from a Pleistocene polar bear genome. Proceedings of the National Academy of Sciences 2022, 119:e2200016119.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J Roy Stat Soc: Ser B (Methodol). 1995;57:289–300.
Salojärvi J, Smolander O-P, Nieminen K, Rajaraman S, Safronov O, Safdari P, Lamminmäki A, Immanen J, Lan T, Tanskanen J. Genome sequencing and population genomic analyses provide insights into the adaptive landscape of silver birch. Nat Genet. 2017;49:904–12.
Chang CC, Chow CC, Tellier LC, Vattikuti S, Purcell SM, Lee JJ. Second-generation PLINK: rising to the challenge of larger and richer datasets. Gigascience 2015, 4:s13742-13015-10047-13748.
Soraggi S, Wiuf C, Albrechtsen A. Powerful inference with the D-statistic on low-coverage whole-genome data. G3: Genes, Genomes, Genetics 2018, 8:551–566.
Klopfenstein DV, Zhang L, Pedersen BS, Ramírez F, Warwick Vesztrocy A, Naldi A, Mungall CJ, Yunes JM, Botvinnik O, Weigel M, Dampier W, Dessimoz C, Flick P, Tang H. GOATOOLS: A Python library for Gene Ontology analyses. Sci Rep. 2018;8:10872.
Jakobsson M, Mayer LA, Bringensparr C, Castro CF, Mohammad R, Johnson P, Ketter T, Accettella D, Amblas D, An L, Arndt JE, Canals M, Casamor JL, Chauché N, Coakley B, Danielson S, Demarte M, Dickson M-L, Dorschel B, Dowdeswell JA, Dreutter S, Fremand AC, Gallant D, Hall JK, Hehemann L, Hodnesdal H, Hong J, Ivaldi R, Kane E, Klaucke I, Krawczyk DW, Kristoffersen Y, Kuipers BR, Millan R, Masetti G, Morlighem M, Noormets R, Prescott MM, Rebesco M, Rignot E, Semiletov I, Tate AJ, Travaglini P, Velicogna I, Weatherall P, Weinrebe W, Willis JK, Wood M, Zarayskaya Y, Zhang T, Zimmermann M, Zinglersen KB. The International Bathymetric Chart of the Arctic Ocean Version 4.0. Sci Data. 2020;7:176.

No competing interests reported.

Download PDF

Editor assigned by journal
11 Sep, 2024
Submission checks completed at journal
05 Sep, 2024
First submitted to journal
04 Sep, 2024

You are reading this latest preprint version

Ancient genomes of Sitka black-tailed deer show evidence for postglacial stepping-stone dispersal along the Pacific Northwest Coast of North America

Status:

Version 1

Abstract

Background:

Results:

Conclusion:

Figures

Introduction

Results

Two distinct maternal lineages inhabited Southeast Alaska during the Holocene

Signatures of mitochondrial-nuclear discordance

Results from exploring signatures of gene flow in further detail

Levels of genetic diversity

FST outliers

Discussion

Deer in British Columbia

The arrival of deer in Southeast Alaska

A postglacial northwards stepping-stone dispersal of black-tailed deer

Conclusion

Material and Methods

Sampling and radiocarbon dating

DNA extraction, PCR amplification, mitochondrial enrichment, and shotgun sequencing

Genome mapping assembly and DNA degradation assessment

Mitochondrial genome analyses

Whole genome analyses

Analyses of population structure and phylogenetic reconstruction

Tests of admixture among Odocoileus hemionus groups

Genetic diversity

FST outliers

Deglaciation and exposure of continental shelves along the Northwest Pacific Coast

Declarations

Conflict of interest statement

Supplementary Tables

Table S3

Author Contribution

Acknowledgment

Data Availability

References

Additional Declarations

Supplementary Files

Status:

Version 1

F_ST outliers

F_ST outliers