Extrinsically reinforced hybrid speciation within Holarctic ermine (Mustela spp.) produces an insular endemic

Refugial isolation during glaciation is an established driver of speciation; however, the opposing role of interglacial population expansion, secondary contact, and gene flow on the diversification process remains less understood. The consequences of glacial cycling on diversity are complex and especially so for archipelago species, which experience dramatic fluctuations in connectivity in response to both lower sea levels during glacial events and increased fragmentation during glacial recession. We test whether extended refugial isolation has led to the divergence of genetically and morphologically distinct species within Holarctic ermine (Mustela erminea), a small cosmopolitan carnivore species that harbours 34 extant subspecies, 14 of which are insular endemics.


| INTRODUC TI ON
Extended isolation in glacial refugia is a documented driver of diversification worldwide (Bennett & Provan, 2008;Hewitt, 1996;Stewart et al., 2010). The cyclic expansion and contraction of glaciers through the Quaternary (2.6 Mya-present) is credited with structuring flora and fauna at the high latitudes (Pielou, 2008), whereby episodes of allopatric refugial isolation lead to genetic and morphological differentiation over time (e.g. mammals, Heaton et al., 1996, Sawyer et al., 2019plants, DeChaine et al., 2014;insects, Ujvárosi et al., 2010). Subsequent periods of glacial recession, however, are often associated with secondary contact between expanding refugial populations. Secondary contact has been primarily viewed as a homogenizing force, with gene flow often erasing differences accumulated between refugial taxa (Gilman & Behm, 2011). Interestingly, established biotic responses to glaciation by terrestrial species-isolation during glacial advances and expansion following deglaciationmay be reversed in island or archipelago systems where isolation and gene flow are further complicated by fluctuations in sea levels which change interisland and island-mainland connectivity (Esselstyn et al., 2009;Fernández-Palacios et al., 2016;Grant & Grant, 2016;Heaney et al., 2005;Sato, 2016). Glaciation typically fractures the ranges of continental species, however, the retention of water in glaciers also lowers sea levels, thereby increasing connectivity between proximal islands through the exposure of shallow continental shelf (Sato, 2016). In turn, glacial recession can release mainland species from refugial isolation, but correspondingly fragment island systems by inundating shallow areas of continental shelf with water, thus creating dispersal barriers and potentially promoting allopatric divergence (Fernández-Palacios et al., 2016). Islands separated by deeper oceanic barriers (e.g. Ireland relative to the other British Isles) or greater geographic distances (e.g. oceanic archipelagos) likely experience fewer bouts of gene flow due to infrequent connectivity with or colonization by mainland source populations.
Climate-associated changes in sea level influence the evolutionary histories of species on both oceanic (Grant & Grant, 2016;Rocha et al., 2016) and near-shore, continental archipelagos (Colella et al., 2018;Sato, 2016;Slager et al., 2020;Whittaker & Fernandez-Palacios, 2007) spanning climate, elevational, and latitudinal gradients (tropical, Caujapé-Castells et al., 2017;temperate, Craw et al., 2017;arctic, high-elevation sky-islands, Folk et al., 2018;Flantua & Hooghiemstra, 2018). Due to the extended distance of oceanic archipelagos from continental source populations, hybridization during low sea levels primarily occurs between established, actively diverging sister taxa that inhabit geographically neighbouring islands that are more fully disconnected at higher sea levels, as seen on the islands of Galapagos and Seychelles (Grant & Grant, 2016;Rocha et al., 2016). In contrast, the flora and fauna of near-shore continental archipelagos experience both periodic connectivity with neighbouring islands and frequent colonization or invasion by mainland dispersers, potentially leading to a more complex history of isolation and contact with multiple different populations, that is less dependent on long-distance dispersal. Higher latitude continental (as opposed to oceanic) archipelago flora and fauna, including that of the Alexander Archipelago of Alaska (USA), Haida Gwaii and Arctic archipelagos of Canada, the Japanese Archipelago, and the British Isles, may be particularly impacted by glacial cycling, as these archipelagos also served as glacial refugia for displaced continental species, and have dynamic histories of connectivity and isolation. The periodic influx of novel genetic material into continental archipelago species through introgressive hybridization with either mainland dispersers or neighbouring island populations may contribute to elevated levels of endemism within continental archipelagos (Grant & Grant, 2014;Larsen et al., 2010).
Patterns of gene flow can now be explored in detail using whole-genome sequencing (WGS), which enables fine-scale detection of historical and contemporary introgression through the reconstruction of demographic patterns in response to glaciation.
Consistent with the coastal refugium hypothesis (Heusser, 1989), classic morphometrics (Eger, 1990) and WGS of the geographically widespread Holarctic ermine (Mustela erminea, King, 1983) have revealed four genetically distinct clades, geographically coincident with three glacial refugia in North America and a fourth that spanned Beringia through the Last Glacial Maximum (LGM, 26.5-19 kya;Colella et al., 2018;Dawson et al., 2014;Fleming & Cook, 2002). One clade, now endemic to two archipelagos along North America's North Pacific Coast (NPC), may represent a distinct species of hybrid origin (Colella et al., 2018). In this case, ancient admixture is hypothesized to have occurred along the NPC during a previous interglacial (~377 Kya; Colella et al., 2018), during secondary contact between the Old World Beringian and New World "East" clades (Fleming & Cook, 2002). Subsequent refugial isolation of the admixed population, potentially within the NPC archipelagos and neighbouring coastline (Colella et al., 2018;Dawson et al., 2014;Eger, 1990), at least since the LGM (~21 Kya) is hypothesized to have promoted allopatric divergence and potentially speciation of Pacific Coast. Significant environmental modification (e.g. industrial-scale logging, mining) has been proposed for a number of these islands, which may elevate the risk of extinction of insular palaeoendemics.

K E Y W O R D S
carnivore, geometric morphometrics, hybrid speciation, reinforcement, speciation with gene flow, stoat this insular population, which may now represent an evolutionary novelty. Following the LGM, sea levels rose in response to glacial recession, and instead of being released from refugia, many NPC populations remained isolated on islands leading to the evolution of substantial intraspecific diversification. Various numbers of ermine subspecies, from 3 to 34, have been recognized (Corbet, 1978;Ellerman & Morrison-Scott, 1951;King, 1983;Hall, 1951), but with no consensus. Given the shared phylogeographic history of numerous NPC endemics (Sawyer et al., 2019), hybridization and isolation (both refugial and insular) acting as an engine for accumulating endemicity is likely a shared process among similarly distributed taxa and among high-latitude archipelagos (Colella et al., 2018;Patiño et al., 2017;Sawyer et al., 2019;Slager et al., 2020).
Whole-genome sequencing remains expensive, particularly for organisms with large genomes (e.g. plants, mammals), and previous sampling limitations (Colella et al., 2018) did not provide for spatially extensive tests of cryptic speciation within the M. erminea complex.
An incomplete understanding of species diversity hinders effective management of insular endemics, including seven ermine subspecies endemic to the NPC (Hall, 1951). Proposed habitat modification, including industrial-scale timber harvests and mining in the Tongass National Forest (USDA, 2007a(USDA, , 2007bStewart, 2016), increases the exigency of characterizing cryptic diversity and novel evolutionary processes within the biogeographically distinct "Sitkan District" (Cook & MacDonald, 2001;Dawson et al., 2007). Reduced-representation genetic and geometric morphometric approaches to species delimitation afford expanded geographic sampling for increased resolution into cryptic species diversity across the complex landscape of the North Pacific. Morphological and genetic changes do not necessarily occur in a regular order (de Quieroz, 1998). Therefore, the inclusion of morphological perspectives in species delimitation enables a more powerful evaluation of congruence across multiple methods and datasets to limit methodological biases (Carstens et al., 2013).
Integrative approaches to species delimitation can help identify distinct population segments (DPS) relevant under the U.S. Endangered Species Act (1973), particularly in organisms that span broad, spatially heterogeneous environments where we might expect profound differences in size due to latitudinal clines (Bergmann, 1847) or insularization (Foster, 1964). The rate of insular morphological and molecular evolution in mammals varies substantially (Barton, 1996;Millien, 2006;Raia & Meiri, 2011), presumably due to the combined effects of isolation and limited effective population sizes, which increase the strength of genetic drift, can alter fitness landscapes and can, ultimately, lead to the evolution of novel phenotypes.
Here, we use the geographically widespread ermine as a model to understand the nuanced role of climate cycling in the accumulation of diversity within continental archipelagos. We integrate geometric morphometric and genetic datasets with expanded sampling from across the Holarctic range of ermine (sampling 27 of 34 nominal subspecies) to test for cryptic diversity. Specifically, we aim to F I G U R E 1 Map of ermine subspecies distributions based on Hall (1951) and IUCN distribution data (www.iucn. org) with individual sampling localities indicated by black dots. (a) Polar view of the northern hemisphere; (b) North America; (c) North America's North Pacific Coast (NPC), geographic limits correspond to the extent-indicator box in (a). The NPC Island clade is distributed across two islands in the Haida Gwaii Archipelago of British Columbia (Graham and Moresby islands; haidarum), on Prince of Wales Island (celenda) in the Alexander Archipelago of Southeast Alaska, and possibly on neighboring Suemez Island (seclusa) determine the species status of the NPC Island clade, which may represent an evolutionary novelty in mammals and provide rare insight into a constructive role for interglacial hybridization in the diversification of archipelago systems that may be broadly relevant to high-latitude refugial archipelagos. all four genetic clades, although sampling is biased towards New World localities. Although the number of ermine subspecies remains contentious, we chose to assess the most diverse estimates (Hall, 1951), hypothesizing that a subset therein (as proposed by Corbet, 1978or King, 1983 would be identified as distinct. Ermine are strongly sexually dimorphic (Hall, 1951;King, 1983); therefore, following Abramov and Baryshnikov (2000) we restricted morphological sampling to males only, as species diagnosis within mustelids remains consistent even when sex is unknown (Bornholdt et al., 2013). Crania were loaned from the University of New Mexico  (Dayan & Simberloff, 1994;Hall, 1951;Peters et al., 2010). Images were collected using a macro lens on a Canon DSLR D5200 fixed 33 cm above the specimen.

| Geometric morphometric methods
Specimens were aligned by eye, with the hard palate parallel to the imaging plane. Orientation aimed to maximize the number of reproducible, homologous landmarks available within the same plane of focus. Images were taken at every degree of focus and digitally stacked into a single high-resolution image using the Helicon software suite (Remote and Focus; www.helic onsoft.com). Homologous ventral landmarks (N = 24), selected based on comparable studies of cranial morphology in carnivores (Goswami et al., 2011;Peters et al., 2010), were digitized in TPSdig (Rohlf, 2006; Supporting Information S1; Figure 2) by the same observer (LMF) to minimize placement error. We estimated landmark variance (morphol.disparity function in the Geomorph package in R v. 3.6.1; Adams et al., 2017;R Core Team, 2017), after landmarking the same specimen image five times on different days to assess placement accuracy. Only the right side of each skull was landmarked to prevent confounding issues associated with fluctuating asymmetry. Unless noted, all statistical analyses were conducted in Geomorph (Adams et al., 2017).
At least one landmark was missing from 106 specimens. For each specimen with missing data, individuals lacked an average of 2.03 landmarks (range 1-7). Missing landmarks were estimated using the estimate.missing function within geographic populations (e.g. island) or subspecies, if geography was not obviously disjointed. A general Procrustes analysis (GPA) eliminated variance caused by rotation, translation, and scale based on the unit centroid size (CS, the square root of the sum of squared distances of all landmarks from their centroid). Nine outliers were identified (plotOutliers), based on distance from the mean shape (above the uppermost and below the lowermost quartiles) and removed from downstream analyses. An additional GPA was run on the final curated dataset, excluding outliers, agnostic of a priori subspecies identifications.
We ran an allometric ANOVA (procD.allometry) on the GPAadjusted coordinates to test for the influence of size (cranium CS, Klingenberg, 2016) on shape variation in male M. erminea (p < .05).
We extracted residual shape variation and ran a Procrustes ANOVA (procD.lm) to test which factors (genetic clade, subspecies, collection locality, or insular/continental collection locality) most influence shape variation. Pairwise Procrustes ANOVAs between subspecies and clades tested for significant differences between groups (p < .05) with the expectation that traditionally defined morphological subspecies would lack significant differentiation once size was F I G U R E 2 The ventral aspect of a landmarked male ermine specimen from the University of Alaska Museum of the North (UAM:Mamm:93279). Landmark abbreviations correspond to the descriptions in Supporting Information S1 eliminated from our analysis. To assess morphological separation among subspecies and genetic clades, we conducted a principal component analysis (PCA). To simplify PCA visualization, we also averaged principal components within each subspecies and within each genetic clade and plotted groupwise averages. We estimated mean shapes (mshape) for each group (subspecies and clade) and used displacement plots (plotRefToTarget) to visualize differences between groups. Mean CS (±SD) was calculated and contrasted between groups as a proxy for overall size (Bookstein, 1996;Slice et al., 1996). A standard linear regression (lm function in base R) estimated the strength of the relationship (R 2 ) between size and latitude.
CS was plotted against collection latitude to visualize agreement with Bergmann's latitudinal rule (Bergmann, 1847).
In addition to testing for differentiation among predefined subspecies and genetic clades, we also ran an uninformed cluster analysis for all New World specimens, for which we had greater sampling depth and subspecies representation, to determine the best partitioning scheme for our morphological data (kmeans function in the flexclust R package; Leisch, 2006). We used the Rand, Jaccard, and Fowlkes-Mallows index (randIndex function; Hubert & Arabie, 1985;Rand, 1971) to compare agreement among kmeans clustering schemes and genetic clusters identified using genetic and WGS data (Colella et al., 2018;Dawson et al., 2014). The Rand Index varies between zero and one, with one indicating that two clustering schemes are identical.
To place New World morphological variation into a broader context, we added limited geometric morphometric coordinate data from 44 Old World ermine, representing eight of 14 Old World subspecies. Differences were visualized using PCA and linear regression to assess the relationship between CS and collection latitude. We contrasted allometric trends between Old and New World ermine and tested for differences in size and shape between subspecies and genetic clades (Old World samples pooled with the Beringia clade) using pairwise ANOVAs, corrected for multiple comparisons.

| Molecular methods
We conducted five runs for each guide tree with different random seeds. Following Berriman et al. (2018), the MCMC chain was run for 2.5 × 10 6 generations, with a 2 × 10 5 burn-in, sampling every 250 generations. Pooled across runs, a posterior probability of > 0.90 was considered significant support for a descendant species.
As an additional delimitation test, we used the bPTP Web server available at: species.h-its.org/ptp/ (Zhang et al., 2013). bPTP generates ML and Bayesian support values following Poisson tree processes. bPTP was run for the ML mitochondrial and nuclear phylogenies generated by Colella et al. (2018) based on 10 ermine WGS.
However, because bPTP is not designed to accommodate multi-locus or SNP data, we then ran each locus generated as part of this study (mitochondrial genomes, ASIP, FES, GHR, HTR1B) independently and as a concatenated supermatrix with corresponding sequences from other Mustela species, as available through GenBank (Supporting Methods). Each bPTP iteration was run against the species-level ML phylogeny generated in RAxML (10,000 iterations, random seed; Stamatakis, 2014) for 100,000 generation with 0.1 burn-in, 100 thinning and a random seed, and convergence was assessed through visual inspection of trace output.

| RE SULTS
Procrustes variance around landmarks averaged 0.0003, indicating high precision in landmark digitization. Ten percent of the total shape variation in New World ermine was correlated with size (Supporting Information S3). After removing allometric variation, a priori subspecies assignment was the most influential variable examined, accounting for almost 16% of size-free variation, followed by collection locality and genetic clade assignment (7% each; Supporting Information S4). Whether a specimen was collected on an island versus from a mainland locality did not significantly influence sizefree shape variation (Supporting Information S4). Although most factors examined contributed significantly to shape variation, 69% of shape variation remained unexplained. Pairwise ANOVAs between the 19 sampled New World subspecies found significant shape and size differentiation (Supporting Information S5). Intraspecific systematic results and discussion are available in the Supporting Information. We found significant morphological differentiation in both size and shape among the four genetic clades (Table 1), with the exception of the East and West clades. Comparison between mean clade shapes (Figures 3 and 4)  Consistent with our ANOVA results, the mean shape of the East and West clades did not differ dramatically (Figure 3). While the palaeoendemic NPC Island clade is the most morphologically distinct (Figures 3 and 4), its phenotype is also most similar to its sister, the Beringia clade (Colella et al., 2018). Both clades share a posteriorly shifted rostral and palatal region and anteroposteriorly compressed braincase, relative to the East clade ( Figure 3). The West clade is morphologically similar to its geographic neighbour and sister, the East clade. Within New World ermine, the Beringia clade has the largest CS (Table 2; Figures 3 and 4), followed by the NPC Island clade, and the West clade has the smallest CS.
Uninformed cluster analyses aligned more closely with genetic clade assignments than subspecific assignments, although neither and M. e. aestiva the smallest sampled, although additional sampling from these regions is necessary to refine these analyses.
After removing shape variation correlated with allometry, we found significant size and shape differentiation (p < .01) between New and Old World ermine (Supporting Information S8-S13). New World ermine have a more elongated braincase, compressed rostrum and anteriorly shifted facial structure relative to the few sampled Old World ermine (Supporting Information S9). For this combined dataset, subspecies assignment had the highest correlation with shape variation of any factor examined (19%), followed by collection locality (9%; Supporting Information S10). Again, 65% of shape

F I G U R E 4 (a) Estimated geographic distributions of
Mustela erminea genetic clades based on the IUCN distribution of the species complex and the distribution of mitochondrial haplotypes (Dawson et al. 2014). Principal component analysis (PCA) of residual shape (or "size free") variation with principal components (PC) one and two averaged by subspecies (b) and genetic clades (c). Colors correspond to genetic clade assignment and subspecies name abbreviations are listed in Table 2 TA B L E 2 Geometric morphometric (GMM) and molecular (Molec.) subspecific (Subspp.) sampling (N), type localities, average group centroid size (CS), ± a standard deviation (SD), and genetic clade assignments  there is mixed support for the morphological differentiation of several Old World subspecies, likely a consequence of low sample sizes.
For example, four Old World subspecies (hibernica, stabilis, aestiva and erminea) exhibit significant size and shape differentiation, while 2 (mongolica and ricinae) show limited or inconsistent differentiation, and others differ only in shape, but not size (ferghanae and kaneii; Supporting Information S5); however, this may be an artefact of small samples sizes.

| Species delimitation
Newly generated sequences (57 mitochondrial  Tests of four species, which recognize both the East and West clades, were also supported (>90% PP; Supporting Information S15).
Genetic data alone parsed four species within the M. erminea species complex; however, tests of geometric morphometric data independently failed to parse the East and West clades as distinct species.
Connectivity and hybridization have been historically viewed as homogenizing forces (Cody, 2006;Whittaker & Fernández-Palacios, 2007); however, episodic bouts of gene flow and geographic isolation may synergistically elevate levels of insular endemism. Hybridization upon secondary contact introduces novel genetic material into small, otherwise isolated populations.
Introgressed variants can accumulate cyclically (high allelic turnover) or gradually overtime depending on the duration of isolation and contact. Introgressed variation may occasionally form the genetic framework necessary for rapid evolution of novel phenotypes, occasionally observed in island systems (Alcala & Vuilleumier, 2014;Caujapé-Castells et al., 2017;Larsen et al., 2010). For ermine, climate-mediated gene flow appears to have acted as a biodiversity "pump," periodically introducing novel variation (April et al., 2013;Avise et al., 1998;Colella et al., 2018;Haffer, 1969;Irl et al., 2017) that was slowly shaped through extended allopatry and a combina- Preble, 1898), as originally described, is a distinct species and represents a specialized case of hybrid speciation in mammals driven by extrinsic reproductive barriers (e.g. geographic isolation; Ottenburghs, 2018). Ancient hybridization between two diverging species is hypothesized to have occurred prior to the LGM (> 300 kya; Colella et al., 2018). We hypothesize that genetic mixing along the NPC was followed by extended geographic isolation of the hybrid population, first in a coastal glacial refugium (through the LGM) and later on at least three islands along the NPC, after sea levels rose in the late Pleistocene. Extended refugial and, now, insular isolation have provided sufficient time for speciation to occur. An influx of novel genetic material, followed by extended allopatry, may have contributed to the rapid divergence of insular endemic Haida ermine (Alcala & Vuilleumier, 2014). Divergence date estimates, based on a log-normal ermine fossil calibration (1.8 Mya, King, 1983) and complete mitochondrial genomes, suggest an initial split between Old World (NPC and Beringia) and New World (East and West) clades around 4.5 Mya. Subsequent divergence of Haida ermine (NPC) from the Beringia clade is dated to 2 Mya (Colella et al., 2018), again, sufficient for speciation to have occurred. The early hybrid history of M. haidarum along the NPC may be analogous to the contemporary ermine hybrid zone observed along the northern Alaska-Canada border today (Colella et al., 2018). The ancient hybrid signature evident in Haida ermine today may persist only as a consequence of subsequent and extended allopatric isolation, while the continental hybrid zone is expected to be more prone to break down over time and will ultimately contribute less to diversification processes.
The geologic complexity of archipelagos may disproportionately preserve evidence of historical introgression, which may otherwise be lost through homogenization upon secondary contact (Currat et al., 2008). Integration of genomic evidence of introgression with more spatially extensive sampling and species delimitation methods clarifies of the varied impacts of range contraction and expansion, gene flow and diversification across the high latitudes, a framework that is essential for anticipating the evolutionary impacts of rapid environmental change now underway and projected into the future (Hope et al., 2013). Additional insular sampling and comparative analyses across varied taxa will better illuminate the generality of climate-driven isolation and contact as a mechanism of diversification. America (Martinkova et al., 2007). In contrast, genetic similarity be-  et al., 1999Sato, 2016). Similar distributional and introgression patterns have been observed in hare of the genus Lepus (Kinoshita et al., 2019). Introgression coincident with climatic cycling has been identified among east Asian mountain hares (Lepus timidus) and Manchurian hares (L. mandshuricus); however, similar to ermine, no introgression was found in Japanese hares (L. brachyurus) on the more isolated island of Honshu (Kinoshita et al., 2019). Interestingly, Corbet (1978) historically recognized only three ermine subspecies throughout the Holarctic: one in Ireland, one on Honshu Island, and one widespread throughout Eurasia and North America. Although Corbet (1978) did not sample ermine from the Alexanderor Haida Gwaii archipelagos, his hypothesis recognizes the evolutionary distinctiveness of refugial continental archipelagos and is consistent with our hypothesis that climate-driven cycles of connectivity and isolation may lead to the accumulation of endemicity in these areas.
Identification of a hybrid carnivore species that is endemic to the Alexander and Haida Gwaii archipelagos suggests that the geographic complexity of continental archipelagos may foster homoploid hybrid speciation or speciation with gene flow at an elevated frequency (i.e. hybridization without a change in chromosome number; Grant & Grant, 2014Larsen et al., 2010;Schumer et al., 2014). Insular hybrid divergence of a non-volant mammal, with limited dispersal capacities relative to bats (Larsen et al., 2010) or birds (Grant & Grant, 2014, suggests that homoploid hybrid speciation may be more common in continental archipelagos, relative to oceanic. Secondary contact can occur more frequently within Ottenburghs (2018) distinguishes between type I ("intrinsic") and type II ("extrinsic") hybrid speciation to account for additional processes driving reproductive isolation (e.g. vicariance) whereby Haida ermine represent an example of extrinsically mediated, type II homoploid hybrid speciation. Ultimately, this cyclic, deep-time process has important implications for continued discoveries of cryptic diversity in these regions and the conservation of endemism across complex landscapes.
Under a morphological species concept (Cronquist, 1978), cranial morphometrics support the recognition of three species within M. erminea; however, under a phylogenetic species concept (Baum & Shaw, 1995;Cracraft, 1983;Donoghue, 1985) both WGS (Colella et al., 2018) Eger, 1990;Hall, 1951), morphology alone does not detect significant divergence between East and West clades. That disparity is consistent with previous morphological investigations that identified clinal variation in size among North American ermine, explained by local adaptation to temperature and precipitation rather than geographic distance (Eger, 1990).
Different patterns among molecules and morphology may be unsurprising for mustelids, which are known to exhibit extensive ecomorphological variation (Dumont et al., 2015;Law et al., 2018).
Therefore, under the general lineage concept of species (Florio et al., 2012;De Queiroz, 2007) Ecomorphological plasticity among mustelids and the role of spatio-temporal and geographic variation in shaping morphological variation are reflected in 30 morphologically distinct ermine subspecies also supported by these analyses (Table 3 and Supporting   Information). This result is better aligned with New World estimates of ermine subspecific diversity proposed by Hall (1951) than the fewer subspecies proposed by Corbet (1978) or King (1983 Cabrera, 1913;King, 1983), which spans much of Eurasia and into Alaska, consistent with sequence data from Kurose et al. (2005). We recommend the nomen Mustela richard- If the NPC was completely ice-covered throughout the LGM (Klein, 1965;Lesnek et al. 2018, Lesnek et al. 2020, there would be a low expectation of endemism in this region. In contrast, frequent discoveries of deeply divergent endemics along the NPC (Barry & Tallmon, 2010;Dawson et al., 2007;Sawyer et al., 2019) provide mounting evidence for the existence of one and possibly multiple coastal refugia along the NPC (Ager, 2019;Baichtal & Carlson, 2010;Carrara et al., 2003Carrara et al., , 2007. Our results suggest that conservation efforts focus on insular endemics (Cook et al., 2006;Sikes & Stockbridge, 2013), as their unique molecular and morphological characteristics may rise to the species level with the analysis of additional data streams. Re-evaluation of cryptic diversity along the NPC is warranted and urgent given anthropogenic changes now underway and planned on numerous coastal islands (USDA, 2007a,b;Stewart, 2016). Old-growth timber harvests within the Tongass National Forest and industrial-scale mining have reduced available forest habitat on these islands (USDA, 2007a,b;Stewart, 2016) and may disproportionately impact insular endemics, such as the Haida ermine, and other non-volant species. Expanding human populations and other impacts (e.g. increasing tourism) along the NPC and the geographic proximity of these islands to the mainland also increases risks of anthropogenic wildlife introductions of invasive species or pathogens that may compete with or immunologically decimate naïve endemics (Ritchie & Johnson, 2009). For example, multiple pathogens common to pets, including canine distemper and parvoviruses, have negatively impacted wild mustelids (Keller et al., 2012;Williams et al., 1988), and the devastating reports of recent TA B L E 3 Proposed revised ermine taxonomy recognizing three evolutionarily independent Mustela species coincident with phylogeography and at least 22 subspecies therein. Twelve additional subspecies require further examination: * indicates low sample sizes in this study, ** indicates unexamined transmission of SARS-CoV-2 from humans to mustelids are yet another potential conservation challenge (Franklin & Bevins, 2020;Gryseels et al., 2020;Molenaar et al., 2020). Given our still rudimentary understanding of much of the diversity in large northern archipelagos, threats to insular endemic diversity should be better quantified to preemptively inform management decisions.

ACK N OWLED G EM ENTS
We dedicate this work to the memory of David R. Klein, an inspirational pioneer in island biogeography and wildlife biology in the North, and especially in the Alexander Archipelago. We thank B McLean for guidance and thoughtful discussion, and Landmark coordinate data (TPS files) and associated images (.JPG) are available from the Dryad Digital Repository: http://dx.doi. org/10.5061d ryad.754vq8n.