Whole-genome analysis of Mustela erminea finds that pulsed hybridization impacts evolution at high latitudes

At high latitudes, climatic shifts hypothetically initiate recurrent episodes of divergence by isolating populations in glacial refugia—ice-free regions that enable terrestrial species persistence. Upon glacial recession, populations subsequently expand and often come into contact with other independently diverging populations, resulting in gene flow. To understand how recurrent periods of isolation and contact may have impacted evolution at high latitudes, we investigated introgression dynamics in the stoat (Mustela erminea), a Holarctic mammalian carnivore, using whole-genome sequences. We identify two spatio-temporally distinct episodes of introgression coincident with large-scale climatic shifts: contemporary introgression in a mainland contact zone and ancient contact ~200 km south of the contemporary zone, in the archipelagos along North America’s North Pacific Coast. Repeated episodes of gene flow highlight the central role of cyclic climates in structuring high-latitude diversity, through refugial divergence and introgressive hybridization. When introgression is followed by allopatric isolation (e.g., insularization) it may ultimately expedite divergence.

G lacial refugia are fragmented pockets of unglaciated land where terrestrial species have persisted and diverged through ice ages, and which have significantly structured biological diversity at high latitudes [1][2][3] . Three primary macrorefugial locations are hypothesized to have enabled the persistence of temperate terrestrial North American species through Pleistocene glacial cycles (>24 glacial/interglacial cycles; 2.6 Mya-11.7 kya). These include an expansive Northern Beringian refugium shared with Asia 4 and two Southern refugia, separated into East and West by the Rocky Mountains or Great Plains/ Mississippi drainage 5,6 . A fourth, smaller refugium (or series of refugia), hypothesized to exist on areas of exposed continental shelf during the last glacial maximum near contemporary North Pacific Coastal archipelagos 7-10 , may have provided another isolated sanctuary for long-term coastal persistence of terrestrial species. Although this refugium has received varying degrees of empirical support, its existence is critical to understanding dynamics of intercontinental biotic exchange and it may have played a primary role in the early peopling of the Americas 11 .
Phylogeography of high-latitude Mustela erminea, the stoat or ermine (Family: Mustelidae) delineates four genetically distinct mitochondrial lineages geographically corresponding to these four refugia 12 . Although glacial isolation appears to have impacted the long-term evolution of this mammalian meso-carnivore, the effects of glacial recession and subsequent secondary contact between refugial lineages remain unexplored. Therefore, we investigated the consequences of cyclic climates and pulsed introgression on genomic diversity in this circumboreal species complex. During inter-glacial warming, divergent lineages expanded from refugial centers until they reached environmental or biological barriers 13 . In some cases, admixture occurred upon secondary contact with other lineages. If M. erminea experienced serial bouts of climatemediated introgression, a repeated signature of introgressive hybridization should be evident in the genomes of refugial descendants 14 .
Contact zones are windows into evolutionary processes 15 now being explored on a genomic scale. We define introgressive hybridization as interbreeding and the movement of alleles between two genetically distinct lineages. Although often difficult to parse introgression from incomplete lineage sorting in closely related taxa or lineages 16 , investigations of divergent lineages, such as those within the M. erminea complex, provide insight into the genomic consequences of introgression. For example, genetic swamping or homogenization can prevent divergence between mixing lineages, while Bateson-Dobzhansky-Muller [17][18][19] and mitonuclear 20 incompatibilities can lead to hybrid breakdown, thereby promoting divergence through reinforcement 21,22 . Although our understanding of hybridization emerges largely from research on plants 23,24 , the genomic era has uncovered substantial reticulate evolution (e.g., introgressive hybridization) in a growing number of animal species 25,26 , including humans 27 . As histories of divergence and speciation-with-gene flow continue to be discovered (e.g., 28,29 ), the implications of hybridization for the conservation of evolutionarily distinct units appear profound 30 .
To test whether genomic structure in stoats reflects refugial origins and to explore the consequences of contact among lineages, we generated ten whole-genome sequences of representatives from each refugial clade. With expanded genomic coverage, we test for signatures of refugial divergence and introgression, characterize the timing of introgression events relative to climatic oscillations, and infer the impact of introgressive hybridization on divergence across the complex landscape of northwestern North America.

Results
Refugial origins and phylogeography. To test our hypothesis that the four mitochondrial stoat lineages originated as the result of isolation and divergence in each of the four major North American refugia, we generated whole-genome sequences for ten M. erminea drawn from each refugial lineage and two spatially disjunct contact zones (the interior border of Alaska-Yukon Territory and the NPC; Fig. 1). We contrasted mitochondrial and nuclear phylogenies against geographic clade distributions, estimated pairwise diversity metrics (pairwise mitochondrial and nuclear divergence, F ST 31 , relatedness 32 ), and used ADMIX-TURE 33 analyses to define populations.
Our results showed strong mitochondrial ( Fig. 1 Tables 1 and 2). We detected substantial genetic divergence between mainland stoat lineages (e.g., 21% average pairwise difference across >3 million autosomal SNPs between Beringia and East clades), comparable to estimates between Beringia stoats and the outgroup species, M. putorius (domestic ferret, 20%), calling into question whether this widespread mustelid is a single species (Supplementary Table 1). Centrality of the Vermont sample is consistent with widespread expansion of the East lineage from a single eastern refugium followed by intra-clade divergence (Fig. 2b). Similar refugial signatures have been detected in other North American species (e.g., martens 34 and black bears 35 ). Additional structure within the East clade (Fig. 2b) may reflect microrefugial diversification 36 . Geographic proximity (Fig. 1) is not predictive of genetic similarity ( Fig. 2b; e.g., REV and BC), consistent with patterns also uncovered in primary stoat prey species such as voles 37 and lemmings 38 .
The sister relationship between NPC and Beringia stoats, combined with substantial genetic divergence (nuclear >10%; mitochondrial >4%; Supplementary Table 1) between these clades and their current geographic distributions, are consistent with the Coastal Refugium Hypothesis 7-10 . Ancestry estimation 33 identified autonomous East and West populations but a single combined Beringia/NPC population (K = 3; Fig. 2c), again highlighting the paleoendemic ancestry of the NPC Island lineage. Pairwise mitochondrial divergence is also particularly elevated for the NPC Island lineage, bordering the species-level divergence threshold suggested for distinct mammal species 39 .
Detecting introgression and determining source populations. Evidence of recurrent introgressive hybridization, whether as a single or recurrent event, should be evident in the genomes of refugial descendants. We hypothesized that signatures of introgression should be present in individuals collected within the Alaska-Yukon Border hybrid zone and expected that stoats from Beringia and East clades are the most likely source populations based on contemporary clade distributions. In addition to contrasting mitochondrial and nuclear phylogenies for discordance, we performed a series of tests for introgression including ADMIXTURE 33 and f3-statistics (Supplementary Table 3) 40,41 to assess shared genetic drift among all possible combinations of two source populations 42 .
We identified mitonuclear discordance for two individuals collected near the Alaska-Yukon border (herein SYT and YTAK, representing Southern and Northern localities along the contact zone, respectively), which can indicate a history of genomic introgression or incomplete lineage sorting. ADMIXTURE identified the same two individuals as admixed between East (K3: SYT 66%, YTAK 74%) and Beringia (K3: SYT 34%, YTAK 19%) clades or as a single, but distinct population, consistent with shared source populations (Fig. 2c). Both Alaska-Yukon samples were positioned between their hypothesized source populations in our PCA results ( Fig. 2a and b), with SYT closer to Beringia and YTAK closer to the East clade, further supporting their mixed ancestry and suggested asymmetric backcrossing distributions. Lastly, f2 statistics, a measure of distance between populations, were also lowest between the putative Alaska-Yukon hybrids and their suspected source populations, further supporting their derivation from the same source populations (Supplementary Table 4).   Fig. 2a and b). These island stoats have been isolated off the NPC (Prince of Wales Island and the neighboring Haida Gwaii archipelago of British Columbia, Canada [not sequenced here]) for at least 10,000 years 12 ; therefore, we suspect they are derived from an ancient hybridization event that occurred during an earlier interglacial period, prior to the contemporary insularization of this clade. We found no evidence of admixture, however, in the Kupreanof Island sample. Admixture date estimates in years based on an effective population size (Ne) of 375k (ref. 12 ) for the North Pacific Coast island (NPC) and the two Alaska-Yukon hybrids (SYT from Southern Yukon and YTAK from the northern Alaska-Yukon border). Err indicates the standard error and Z indicates the Z-score 8% of the total variation and PC2 11.8%. c ADMIXTURE plots for K2-7. K = 3 best matches our expectation based on the geographic distributions of refugial clades and is consistent with admixed ancestry for SYT and YTAK samples: (1) East samples as a distinct group, (2) a West population, and (3) a combined Beringia-NPC population with two hybrids admixed from these groups. K = 6 has the lowest cross-validation score for our data and shows an identical pattern to K3 but with additional substructure within the East clade, indicative of microrefugial substructure near Appalachia and consistent with PCA results suggesting East expansion from a single refugium Testing for asymmetric introgression. Biased or asymmetric introgression may indicate a selective advantage, the existence of genetic incompatibilities, or a demographic imbalance. Alternatively, equal contributions from both parental source populations may suggest that there are minimal genetic incompatibilities between hybridizing stoat clades, as we would expect within a single species. For the NPC Island clade, biased allelic overlap with one source population may shed light on interglacial colonization dynamics, as the island may have been colonized by both Beringia and East ancestors who subsequently interbred. We used F4-statistics, similar to D-statistics 43,44 and ABBA/BABA 45 , to highlight differential backcrossing histories among lineages and used admixture graph fitting 46 to estimate the topological placement of admixed samples relative to their potential sources. First, our f4 results highlight the influence of repeated glaciation in North America as reflected in multiple episodes of gene flow between lineages originating in independent refugia. Allele sharing (f4) reflected phylogenetic relationships, with stepwise genetic distance from the ancestral Beringia clade increasing with greater geographic distance from Beringia ( Fig. 1; Supplementary Tables 5 and 6), a pattern consistent with more frequent or prolonged contact between geographically proximal refugial lineages.
Admixture graph fitting identified SYT as the result of Beringia and East admixture, while YTAK was produced from an ancestral admixture event between East and a shared ancestor of all stoat lineages ( Fig. 3a and b; Supplementary Table 7). The vast majority of f4 comparisons (Supplementary Tables 5 and 6) were highly significant (Z-score >5), reflecting substantial intra-specific divergence within M. erminea. Biased allelic overlap suggested that YTAK shares more alleles with the East source, whereas SYT is more backcrossed with the Beringia source (   All Z-scores are positive for ease of interpretation and significant (>5). Rotating the order of a pair (e.g., YZ to ZY) switches the sign +/− of both the f4-statistic and the Z-score but does not change the relationships and Beringia lineages may have minimal negative fitness consequences (e.g., reduced fertility, viability) 47 . These preliminary inferences require additional sequencing of hybrids from the interior Alaska-Yukon hybrid zone to characterize the specific genomic regions that may be crossing semi-permeable lineage boundaries and to rigorously test for asymmetric introgression, which could be driven by selection (e.g, genetic incompatibilities, Darwin's corollary to Haldane's rule 48 ), drift, or demographic parameters. If gene flow is neutral, extended periods of contact may promote fusion between East and Beringian lineages. However, the contemporary persistence of both Old World Beringian and New World East lineages in North America suggests that gene flow during brief periods of interglacial contact through the Pleistocene were insufficient to homogenize refugial divergence.
Biased allele sharing between NPC and Beringia ( Table 2) is consistent with phylogenetic relationships (Fig. 1d) and three possible evolutionary histories: (1) differential introgression (selection) favoring Beringia alleles early in the NPC Island clade's origin; (2) extended or more frequent connection between Beringia and NPC refugia throughout the Pleistocene; or (3) stochasticity, potentially compounded by founder effects or bottlenecks.
Timing introgressive events. Given the pattern of North American glaciation 49 , the current geographic isolation of the NPC Island clade, and the limited dispersal abilities of insular M. erminea 50 , we suspect that hybridization is not ongoing on NPC islands. In contrast, based on the current distribution of mitochondrial diversity 12 and our understanding of postglacial population expansion, we expect that stoats from East and Beringia clades are actively interbreeding along the Alaska-Yukon border. To approximate the timing of admixture events, we converted drift unit branch lengths (D), output from MixMapper 51 , to absolute time (years) using the formula D ≈ 1−e −t/2Ne (solved for t generations) 52 .
The estimated time of admixture for Alaska-Yukon hybrids did not differ from zero ( Table 1), suggesting that the Alaska-Yukon Territory border is an active hybrid zone between the leadingedge of Beringia and East stoat clades. In contrast, NPC admixture is estimated at up to 394 kya (Table 1; Supplementary  Table 8), prior to the Wisconsin glacial (30-80kya 53,54 ) and Sangamonian interglacial (80-140kya 53,54 ). The weaker NPC hybrid signal relative to contemporary Alaska-Yukon hybrids (f3 Z-score), the insular isolation of this clade for at least the past 10,000 years 12,55 , and the short generation time of stoats (2 years 56 ) suggest that NPC Island admixture is more ancient than that occurring in the mainland Alaska-Yukon zone. We hypothesize that the hybrid origin of the NPC clade predated a period of glaciation and refugial isolation, but this point requires further exploration with expanded sampling and tests of linkage disequilibrium. If NPC stoats continue to evolve in isolation, our results suggest that continued divergence in allopatry may ultimately lead to speciation.
Based on our results, NPC islands were colonized either by admixed individuals or independent invasions of both parental lineages (East, Beringia) that subsequently hybridized. Introgressive hybridization, preceding divergence in isolation, may have disproportionately elevated the evolutionary distance between NPC Island stoats and other lineages (Figs. 1, 2, and 3d, e; Supplementary Tables 1-4). Isolation of a genetically recombinant population, at least through the Last Glacial Maximum, in a coastal refugium reinforced an evolutionary trajectory independent from either source population (Supplementary Fig. 3). Further, low amplicon variation in NPC stoats 12 suggests that island founder effects or a bottleneck 5 occurred among NPC stoats in the coastal refugium 57 . Drift may impact small populations and may have further accelerated NPC stoat differentiation.
Historical demography. We expected glacial cycling (periods of isolation/contact) to have impacted the historical demography of each refugial stoat clade differentially. Further, we hypothesized that individuals with a history of introgression will exhibit elevated effective population sizes (N e ), with a distribution most similar to that of their most-backcrossed source population (East for YTAK, and Beringia for SYT). Last, we expected cladogenesis to coincide with refugial isolation during glacial periods and demographic troughs, while N e is expected to rise during interglacials as a consequence of population expansion.
Stoat demographic histories (PSMC 58 ; Fig. 3d, e) showed clade-level responses to cyclic glaciation, highlighting the central role of large-scale climatic shifts in shaping diversity at high latitudes through episodic isolation and contact. Overall low, but relatively stable N e in East and West clades reflected similar demographic responses to New World glacial cycling, where the West refugium was likely smaller than the East, as reflected in smaller N e (Fig. 3d). Unlike Nearctic populations, Western Beringia experienced substantial expansion historically, corresponding to reduced glaciation in eastern Eurasia and the larger size of Beringia compared to the NPC refugium. Dates of cladogenesis mirrored estimates from the calibrated mitochondrial genome phylogeny ( Supplementary Fig. 1), with inter-clade divergence beginning around 2 Mya when Beringia split from the other stoat lineages, followed by West (1.3 Mya) divergence, and most recently Island/East (0.9 Mya) cladogenesis (Fig. 3d). However, the hybrid ancestry of the NPC island sample suggest that cladogenesis estimates based on N e may be unreliable. As expected, when hybrid distributions are juxtaposed against their parental populations (Fig. 3e), they mirror that of their mostbackcrossed source population, but with an upward shifted N e 59 . If the increase in hybrid N e occurs at the time of hybridization between parental populations 59 , SYT and YTAK admixture may be older than predicted based on admixture date estimates alone.
Combined with diversity metrics (Supplementary Fig. 2; Supplementary Tables 1-4), highly distinct demographic histories suggest taxonomic revision of M. erminea may be warranted, pending wider geographic sampling. This is a particularly compelling example of cryptic diversification in carnivores when contrasted against weakly differentiated PSMC demographic histories uncovered in some canid 60 (wolf, dingo, domestic dog) and felid 61 (tiger, leopard, cheetah) genomes.

Discussion
We document multiple temporally and spatially offset bouts of introgressive contact between divergent stoat lineages, corresponding to distinct phases of glacial cycling. Therefore, we predict that other higher latitude species with shared biogeographic histories may also exhibit genomic signatures of repeated introgression corresponding to cyclic climates. This process of recurrent introgressive hybridization between geographically proximal lineages, followed by refugial or insular isolation and divergence, may contribute to elevated levels of island endemism, particularly for continental archipelagos adjacent to glacially dynamic mainland areas.
When backcrossing occurs into both parental populations, as in the Alaska-Yukon hybrid zone, introgression may be stochastic or selectively neutral. Consequently, extended periods of contact may lead to homogenization 18,62 and result in a net loss of regional diversity 63 , such that the divergent lineages evident today (East, Beringia) may be ephemeral 64 . The narrow Alaska-Yukon hybrid zone and strong East-Beringia genomic differentiation, however, suggests that such peripheral hybrid populations are ephemeral and minimally contribute to long-term evolution of high latitude diversity. If introgression events are temporary, genomic evidence of serial hybridization offers a model for detecting refugia and tracking colonization routes among high latitude species. Here, we provide yet another example of introgressive hybridization in mammals, highlighting both the prevalence and varied role of introgression in mammalian evolution. Our results demonstrate that, for geographically widespread species, the long-term evolutionary role of introgression often may be inconsequential relative to larger biogeographical processes, such as island vicariance, refugial isolation, or rapid environmental change.
In some situations, however, a population bottleneck in a small refugium or in response to insularization can enable the temporary persistence of less fit recombinant genotypes 65 . Even a temporary reduction in fitness can enable traversal between fitness-landscape optima and promulgate alternative evolutionary trajectories, which can permanently impact the genomic structure of high latitude populations and potentially lead to speciation 65,66 . Genomic evidence of hybridization, followed by divergence in geographic isolation, is consistent with a hybrid origin for the NPC lineage. Though still incipient in this case, the process of introgression followed by allopatry may lead to hybrid speciation 67 . Hybridization can catalyze the evolution of reproductive barriers 68 or expedite divergence through the rapid generation of novel genetic recombinants. Although the mechanisms underlying this special case of divergence with gene flow remain elusive, a plausible precondition for speciation is the separation of lineages into allopatric ranges, which may enhance the genetic divergence of a recombinant population, even if introgression is not directly implicated in speciation. Therefore, NPC stoats potentially highlight the integral role of post-hybridization allopatric divergence in the evolutionary histories of island and refugial taxa, but expanded insular sampling is required to fully test this mechanism of divergence.
It has long been argued that pre-zygotic isolation is reinforced when diverging groups return to parapatry and hybridize 18,20 , although this idea is not without controversy. Under this scenario, hybridization between incipient species can result in unfit hybrids and natural selection might act against hybridization, reinforcing pre-zygotic isolation mechanisms, eventually leading to speciation and avoiding maladaptive hybridization 18,69 . In this fashion, continuing and progressive reinforcement of pre-zygotic mechanisms would ultimately result in speciation and complete isolation between lineages 18,70 , if hybridization and selective pressures on both lineages are symmetrical 71 . With an active and plausibly ongoing hybrid zone, along with more ancient island admixture, stoats provide a unique opportunity to investigate the impact of symmetrical versus asymmetrical introgression using genomics.
With multiple waves of mainland colonizers corresponding to climatic cycles and potentially adding novel genetic elements to insular genomic architectures, the process of repeated, climatemediated hybridization may have played a role in diversification across other land-bridge island systems. Irish stoats, for example, exhibit elevated genetic diversity compared to neighboring populations in Great Britain, despite persisting on a smaller island farther from a continental source 13,72 . In this parallel system, the English Channel separating Britain from the European mainland has been serially exposed through Pleistocene glacial cycles 73 , effectively homogenizing populations across the Channel. On the other hand, the Irish Sea, nearly twice the depth of the English Channel, further separates Ireland from Britain, presenting a formidable barrier to gene flow that allowed only occasional influx of novel alleles from mainland or British immigrants during extreme glacial events. This mechanism of post-hybridization allopatric divergence is not restricted to land bridge or continental archipelagos. Oceanic island populations of Darwin's storied finches on the Galapagos Islands also were impacted by climate-mediated pulsed hybridization that was coincident with changes in sea level due to glaciation 74 .
Admixture in this insular endemic (NPC Island stoats, subspecies: M. e. haidarum) poses a particular challenge for the U.S. Endangered Species Act, which currently denies protection for hybrid taxa 30,75 . Mustela erminea haidarum is federally protected in Canada on the Haida Gwaii Islands of British Columbia 56 , but also occurs on Prince of Wales Island, Alaska, where it has no protection. This study is congruent with previous work noting that this island-restricted stoat subspecies is genetically 12,76 and morphologically 77 distinctive, but levels of management protection at the federal level differ profoundly between Canada and the United States.
Genomic signatures of repeated admixture and refugial isolation among North American stoats highlight the episodic role of Late Quaternary climate cycles in generating and structuring high latitude diversity. We demonstrate two instances of gene flow between stoat lineages coincident with interglacial periods: one contemporaneous along the interior Alaska-Yukon border and a second, more ancient event, occurring along the NPC that left a signature that persisted through multiple glacial cycles. Each hybridization event exemplifies an antagonistic outcome of hybridization (e.g., homogenization versus diversity generation) within a single species complex, highlighting the multifaceted role of introgression in mammalian evolution and the varied temporal and geographic scales at which this process unfolds. At high latitudes, where dramatic range shifts in response to glacial cycling afforded recurring opportunities for contact and isolation, ephemeral hybridization has repeatedly shaped biological diversity 2,78,79 .

Methods
Sequencing, assembly and post-processing. Contact zones and populations were defined based on expanded mitochondrial sampling 12 . Ten whole-genome sequences for representatives from each refugial stoat clade, and with additional sampling in suspected contact zones (Fig. 1 Table 9). Liver subsamples were obtained post-mortem from the University of New Mexico's Museum of Southwestern Biology in compliance with ethical regulations (IACUC #16-200526-MC). We used a DNeasy Blood and Tissue Kit (QIAGEN, USA) extraction. Nextgeneration sequencing libraries were prepared using the Illumina TruSeq DNA Sample Prep Kit. Our genome assembly pipeline followed Lan et al. 79 with reads examined using FastQC 80 and adapter sequences and sex chromosomes removed (Trimmomatic v0.33 (ref. 81 )). Quality scales were compared across samples and rescaled to Illumina 1.8 + phred + 33 quality score, when necessary, to maintain consistency (Picard v1.9: http://picard.sourceforge.net; SAMTools 82 ). Reads were mapped to the domestic ferret genome 83 (Mustela putorius furo) using the Burrows-Wheeler aligner 84 (BWA). An additional BWA iteration extracted mitochondrial genomes. Final depth of coverage ranged from 10 to 61× (Supplementary Table 9). PCR duplicates were removed using the MarkDuplicates tools from the Picard suite. Nuclear and mitochondrial genome consensus sequences were called using mpileup in SAMtools. SNPs were called using HaplotypeCaller in the Genomic Analysis Toolkit 85 (GATK) for all stoats and again against each of two outgroups (M. putorius and Martes americana) and filtered by minimum depth (minDP = 2), genotype quality (minGQ = 30), minimum minor allele frequency (MAF = 0.1), and scaffold size and location (1 Mb), and then private alleles and indels were removed using VCFtools 86 (Supplementary Table 10). We intentionally selected an MAF of 0.1 to remove singletons (e.g., individual-specific, rare mutations), and thereby any potential sequencing errors, which are more common in low-coverage genomes. Although singletons can be informative when analyzing within population variation, our analyses focus on shared alleles (e.g., f-statistics) to identify introgression and backcrossing histories, therefore singletons, since they are not shared, do not provide information about allelic overlap among populations. Our minDP (2) was set as one-third of the coverage of our lowest coverage samples (e.g., roughly one-third of 8×) as recommended by the PSMC manual. Also, 2 is the minimum requirement for detecting heterozygotes. Format conversions (vcf, ped, bed) were conducted using PLINK 87 . Missing data were removed (--max-missing) based on analysis specifications. SNPs were spaced (1 per 50 bp window) to account for linkage disequilibrium (e.g., 25) and sorted into 46 'pseudo-chromosomes' to enable the application of human-specific analyses on a non-model system with 40 chromosomes.
Phylogenetic methods. jModeltest 88,89 estimated the most appropriate model of evolution for mitochondrial genomes and autosomal SNPs, with and without outgroup sequences (M. putorius and Martes americana). We generated phylogenies for the entire mitochondrial genome (calibrated/uncalibrated) and autosomal nuclear SNPs (uncalibrated) using RAxML 90 (GTRCAT model, 10,000 generations, random starting seed) and Beast2 91 (1 million generations, 4 chains, 2 runs, sampling frequency of 1000). It should be noted that phylogeny inference using highly variable data (e.g., SNPs) can induce acquisition bias resulting in longer branch lengths 92 , however, other divergence metrics (e.g., PCA, relatedness, F ST ) provide additional evidence of substantial divergence between refugial stoat clades. The oldest known M. erminea fossil 93 95 . To accommodate genetic exchange between populations of the same species, TreeMix uses genome-wide allele frequency data to infer historical relationships among populations to maximize the composite likelihood of the sample covariance matrix (W hat ), while modeling both cladogenesis and migration. TreeMix iterates over a standard bifurcating topology to identify a series of migration edges that most increase the likelihood score. TreeMix was run without default sample size corrections (-noss), because multiple populations were comprised of a single individual. We limited the number of migration events (-m) to three, to examine only the most influential events among the few groups examined. TreeMix simulations with M. putorius as the outgroup identified ancient geneflow between Beringia stoat ancestors and M. putorius. Simulations were re-run with M. americana as an outgroup to assess whether long-branch attraction between distantly diverged M. erminea and M. americana (>11Mya 96 ) impacts migration inference in TreeMix, but results were consistent, regardless of the outgroup used. After TreeMix identified gene flow between M. erminea and M. putorius ( Supplementary Fig. 4) 33 analysis, sites with >80% missing data were removed and the lowest crossvalidation (cv = 10-and 5-fold) score identified the most appropriate number of populations (K = 6) by iteratively leaving a sample out and reexamining the partitioning of genetic structure among the remaining samples.  Fig. 4). f2 statistics [f 2 (Pop 1 , Pop 2 )] quantify drift through summary statistics (e.g., allele frequencies, heterozygosity, covariance, and the probability of two lineages coalescing), with larger values indicating greater divergence 42 (Supplementary Table 4). f2 statistics were generated from the compute_moment_stats and compute_most_additive_trees functions in MixMapper with 1000 bootstraps and SNP blocks of 50 (1 per 50 bp). f3-statistics explicitly test for admixture (AdmixTools, 3PopTest) and considered all permutations where Source 1 , Source 2 , and the Target samples came from different populations. Populations (e.g., clades) were defined by phylogenetic relationships (Fig. 1) and ADMIXTURE results (Fig. 2c). A positive f3 value does not necessarily indicate the absence of admixture 42 . Small sample sizes and the absence of a linkage map for M. erminea prevents linkage disequilibrium-based estimates of N e and refined dating of admixture events.
To parse the backcrossing history of each hybrid sample, f4-statistics were used. We used block-jackknifing to accommodate non-independence between loci 42 . Although f-statistics alone cannot deduce the direction of gene flow in a system, admixture graph fitting can test whether a proposed evolutionary model fits the data well 45,51 . Admixture Graph 46 fit hybrid individuals into a non-admixed tree topology ( Supplementary Fig. 5a) based on f4 results (Supplementary Tables 5 and  6). We ran an additional Admixture Graph simulation on the optimal (lowest minimal error) 4-taxa non-admixed backbone topology identified by Admixture Graph (Supplementary Fig. 5b; Supplementary Table 12) and again for Alaska-Yukon hybrids onto a 5-taxon backbone ( Supplementary Fig. 5c) to evaluate how inferred relationships changed (Supplementary Fig. 6). We compared individuals from each population against members of alternative populations (e.g., f 4 (Outgroup, NPC; East, Beringia)). We also compared hybrid individuals (as identified by f3-statistics) against individuals from 'pure' populations (e.g., f 4 (Outgroup, Hybrid; East, Beringia)) and admixed populations (e.g., f 4 (East, SYT; Beringia; YTAK)) to decipher the backcrossing histories of hybrid samples and characterize patterns of gene flow across populations. f4statistics are negative (Z-score ≤−5) if there is more allelic overlap between X and Y than between X and Z since the evolutionary split between Y and Z, and positive (Z-score ≥5) if there has been more recent gene flow between X and Z than between X and Y. f4 results are in Table 2, Supplementary Tables 5 and 6, and Supplementary Fig. 7. Dawson et al.'s 12 N e Model 2 using multi-locus Approximate Bayesian Computation approximates mean clade N e at: 250k for Beringia, 50k for NPC, 125k for East, and 75k for West; therefore, the combined N e for NPC, East, and Beringia populations is 375k. D is tightly correlated with N e and additional efforts to refine historical population sizes (Pleistocene, Pliocene) should improve the resolution of our admixture date estimates. Generational output was converted to years using a M. erminea generation time of 2 years 56 .
Historical demography. PSMC 58 was used on consensus genomic sequence data (using both complete and down-sampled sequence data) to characterize historical demography (Fig. 3d, e) by examining heterozygosity densities in a 100 bp sliding windows across the genome. PSMC was run twice for each individual, once utilizing all mapped sequence data (coverage 10-61×) and again on data down-sampled to 10× coverage (matching the lowest coverage sample) using the DownsampleSam tool (Picard). Results were scaled by a general mammalian mutation rate (2.2 × 10 −9 per base pair per year 98 ) and stoat generation time, resulting in distributions of N e through time (Fig. 3d, e). One hundred PSMC bootstrap replicates were performed and plotted for both fullcoverage and down-sampled data to confirm consistent distributional shapes and enable comparison across individuals, as PSMC is sensitive to variation in coverage depth (ideal coverage 18× (ref. 99 )).