North American water shrews, which have traditionally included Sorex alaskanus, S. bendirii, and S. palustris, are widely distributed through Nearctic boreal forests and adapted for life in semiaquatic environments. Molecular mitochondrial signatures for these species have recorded an evolutionary history with variable levels of regional divergence, suggesting a strong role of Quaternary environmental change in speciation processes. We expanded molecular analyses, including more-comprehensive rangewide sampling of specimens representing North American water shrew taxa, except S. alaskanus, and sequencing of 4 independent loci from the nuclear and mitochondrial genomes. We investigated relative divergence of insular populations along the North Pacific Coast, and newly recognized diversity from southwestern montane locations, potentially representing refugial isolates. Congruent independent genealogies, lack of definitive evidence for contemporary gene flow, and high support from coalescent species trees indicated differentiation of 4 major geographic lineages over multiple glacial cycles of the late Quaternary, similar to a growing number of boreal taxa. Limited divergence of insular populations suggested colonization following the last glacial. Characterization of southwestern montane diversity will require further sampling but divergence over multiple loci is indicative of a relictual sky-island fauna. We have reviewed and revised North American water shrew taxonomy including the recognition of 3 species within what was previously known as S. palustris. The possibility of gene flow between most distantly related North American water shrew lineages coupled with unresolved early diversification of this group and other sibling species reflects a complex but potentially productive system for investigating speciation processes.
Considering an accelerating loss of biodiversity, expanding human populations, and warming climate trends, understanding how species have responded to changing environments in the past can inform future conservation and management (Parmesan 2006). Many extant species have evolved as a consequence of repeated climate fluctuations through the late Quaternary (Hewitt 1996). With the advent of multilocus molecular techniques, increasingly accurate estimates of evolutionary rates, and extensive archives of contemporary and fossil specimens, evidence is accumulating for common timing of diversification and shared spatial extent of lineages over a broad range of taxa (Brito and Edwards 2009; Ho et al. 2011; Hope et al. 2013a). Terrestrial species responded to change to varying degrees by either persisting in a given area or shifting to follow preferred conditions (Riddle 1996). Although species at lower latitude may have experienced fragmentation and local population expansion–contraction through climate cycles, they often did not drastically shift their range. Desert faunas of the North American Southwest therefore retained a genetic legacy that indicates deep divergence in response to geological events of the Miocene and Pliocene (15–2.5 million years ago [mya]—Riddle 1996; Hafner and Riddle 2011). In contrast, temperate and high-latitude species usually experienced (often obligatory) range shifts that resulted in relatively shallow gene coalescence due to recent speciation or, perhaps more commonly, a loss of historic diversity due to both extirpated populations and pioneer expansion (Hewitt 1996). Discrepancy between deeper fossil histories and shallower genetic structure attests to the loss of some of the genetic record in higher-latitude systems (Dawson et al. 2014).
Shrews of the genus Sorex are a good example of a latitudinal gradient in genetic diversity that reflects deeper evolutionary histories in low-latitude species (Esteva et al. 2010) and shallower genetic differentiation within lineages at higher latitudes (Hope et al. 2013a). This gradient also occurs within species or complexes of closely related species that span broad latitudinal ranges (e.g., Demboski and Cook 2001; Hope et al. 2012). Because much of North America periodically was covered by continental ice sheets during glacial phases, most northern terrestrial species recolonized previously glaciated areas during each subsequent interglacial leading to genetically depauperate northern populations following pioneer expansion. In contrast, populations south of ice sheets retained higher genetic diversity (Lessa et al. 2003), a process accentuated by repeated fragmentation into isolated sky-island populations on southern mountaintops (Knowles 2001; Galbreath et al. 2009, 2010). Such processes exemplify diversification of boreal species associated with temperate forest biomes (Arbogast and Kenagy 2001), the predominant habitats occupied by Sorex. However, there are exceptions to these general patterns, most often realized as northern tundra-associated species (Fedorov et al. 2003; Hope et al. 2013b) or refugial populations confined at high latitudes to discrete periglacial regions (Cook et al. 2001). Refugial populations often are identified by their distinct genetic signatures.
The North American water shrew complex consists of 3 species with a combined distribution that approximates the extent of boreal forest in North America (Fig. 1; Hall 1981). Sorex bendirii (marsh shrew; hereafter Coastal group) is endemic to the western coastal lowlands and Olympic and Cascade mountains of the westernmost United States extending as far north as southwestern British Columbia, Canada (Pattie 1973). O'Neill et al. (2005) found no evidence for significant genetic divergence among 2 of the 3 morphological subspecies of S. bendirii, based on mitochondrial DNA (mtDNA). S. palustris (American water shrew) occurs from central Alaska south to eastern Arizona and northern New Mexico and from Nova Scotia, Canada, in the east to a parapatric distribution with S. bendirii in the west (Fig. 1; Beneski and Stinson 1987). The 3rd species, S. alaskanus (Glacier Bay water shrew), is known from 3 specimens from Glacier Bay, Alaska, and is not included in the current study; the taxonomic validity of this species is suspect (Denver Museum of Nature & Science, K. Hildebrandt, pers. comm.; Hutterer 2005; MacDonald and Cook 2009).
Fig. 1.
Map of the study area indicating general specimen localities for samples of North American water shrews included in the present study. Locality numbers correspond to those listed in Appendix I. Left inset provides finer detail of western localities. Right inset shows polygons representing the geographic sampling of major groups for the current study, fossil localities (black dots; adapted from Himes and Kenagy 2010), and question marks indicating areas where further sampling is required to clarify lineage distributions.
All North American water shrews belong within the “Sorex vagrans complex” of shrews (hereafter vagrans complex), a diverse species complex with unresolved systematic relationships (Findley 1955; Carraway 1990; Demboski and Cook 2001; Shafer and Stewart 2007). Within S. palustris, recent studies have reported substantial mtDNA variation (O'Neill et al. 2005; Himes and Kenagy 2010; Mycroft et al. 2011). Four geographic groups largely are congruent with at least 6 of the 9 described subspecies. Vancouver Island, Canada, supports an isolated subspecies, S. p. brooksi (hereafter V.I. group), which was found to be minimally divergent from the widespread S. p. navigator (hereafter Cordilleran group), a subspecies occupying western North America and spanning the latitudinal extent of the species (O'Neill et al. 2005). Central Canada and the northern midwestern states of the United States support S. p. palustris (hereafter Boreal group). Cordilleran and Boreal groups may be parapatric through central Alberta, Canada, and northeastern British Columbia (Fig. 1), and are polyphyletic with respect to other species of the vagrans complex (O'Neill et al. 2005; Mycroft et al. 2011). Finally, S. p. albibarbis and S. p. gloveralleni (hereafter collectively Eastern group, including S. p. labradorensis) have similar mtDNA, but together were genetically distinct from the Boreal group. Eastern and Boreal groups exhibited moderate mtDNA sequence divergence (Mycroft et al. 2011). All species of North American water shrew are semiaquatic; however, S. bendirii occupies a range of marsh, wetland, and riparian habitats, whereas S. palustris is more strongly associated with forest streams, lakes, and riparian corridors, and rarely found far from clear running or open water (van Zyll de Jong 1983).
North American water shrews represent an excellent group for phylogeographic investigation through the boreal Nearctic. Fossil localities for these species indicate a wide historical range, including areas not currently occupied (Kurtén and Anderson 1980; Fig. 1). Southern populations now occupy isolated sky islands (Smith 1993), whereas coastal island populations along the North Pacific Coast may represent refugial lineages as discovered in other species, or reflect episodes of relatively recent colonization from mainland populations (Nagorsen 1996; Nagorsen and Keddie 2000; Brown and Hebda 2003). Most of the current range of North American water shrews was covered by continental ice sheets through the last glacial phase until following the Last Glacial Maximum (∼21 thousand years ago [kya]), indicating extensive range expansion, although sources for recolonization remain unresolved.
This phylogeographic assessment addresses questions related to spatial and temporal evolutionary processes in North America. Our analyses complement previous assessments by developing a rangewide, multilocus perspective to review systematic relationships and relative divergence within the North American water shrew complex using a Bayesian coalescent framework for species-tree estimation; address the timing of isolation of insular populations along the North Pacific Coast in relation to a glacial climate regime and potential extended refugial isolation; contrast variable demographic responses to recent environmental processes both within and among major North American water shrew lineages; describe novel genetic diversity among North American water shrews that may inform future management of relictual faunas; and report possible contemporary gene flow between major North American water shrew lineages and discuss future research potential of hybrid dynamics.
Materials and Methods
Sampling and sequencing
We used 205 specimens (Appendix I) including 14 S. bendirii (12 S. b. bendirii and 2 S. b. palmeri), 175 S. palustris (3 S. p. albibarbis, 25 S. p. brooksi, 5 S. p. gloveralleni, 1 S. p. labradorensis, 125 S. p. navigator, and 16 S. p. palustris), and 10 outgroup taxa (1 S. cinereus, 3 S. fumeus, 1 S. haydeni, 5 S. monticolus, 1 S. ornatus, 1 S. pacificus, 1 S. rohweri, 1 S. sonomae, 1 S. trowbridgii, and 1 S. vagrans). North American water shrew specimens represented 73 general collecting localities (Fig. 1; 140 specific localities; Appendix I). For 118 of these specimens, sequences for at least 1 gene were retrieved from GenBank. We sequenced the remaining 87 specimens for 1–4 independent loci including 639–1,140 base pairs (bp) of mtDNA cytochrome-b gene (Cytb; n = 201), 434 bp of nuclear apolipoprotein B exon (ApoB; n = 89), 282 bp of interferon 1 intron (IFN1; n = 88), and 286 bp of myosin heavy-chain 2 intron (MYH2; n = 76).
We extracted DNA from frozen heart, kidney, spleen, or muscle tissue (stored at −80°C or ethanol preserved). Genomic DNA was obtained through standard salt extraction with sequences obtained through polymerase chain reaction and cycle sequencing methods. Primers for amplification of each locus included MSB05/MSB14 (Cytb—Hope et al. 2010), ApoBF/ApoBR (ApoB—Dubey et al. 2007), and IFN1F/IFN1R and MYH2F/MYH2R (IFN1 and MYH2, respectively—Lyons et al. 1997).
For all loci, polymerase chain reaction reagents and conditions were: 1 μl of DNA template (variable concentration), 1.5 μL each of deoxynucleoside triphosphates (10 mM), MgCl (25 mM), 10x polymerase chain reaction buffer, bovine serum albumin (1%), 0.2 μl of each primer (10 mM), 0.08 μl of AmpiTaq DNA polymerase (Applied Biosystems, Foster City, California), and 7.52 μL of double-distilled H2O to total 15-μL reactions. Polymerase chain reaction included initial denaturation at 94°C for 6 min, followed by 40 cycles of denaturation at 94°C for 25 s, annealing at 50°C for 30 s, extension at 72°C for 1 min, and final extension at 72°C for 5 min, with cooling at 15°C for 10 min. Polymerase chain reaction products were cleaned using ExoSAP-IT (Affymetrix, Santa Clara, California) diluted 1:10. Reagents for each 12-μl cycle sequencing reaction included 2.5 μl of H2O, 2 μl of BigDye buffer (Applied Biosystems), 1 μl of primer (10 mM), 0.5 μl of ABI BigDye version 3.1 (Applied Biosystems), and 5 μl of DNA (variable concentration). We followed cycle sequencing protocols of Platt et al. (2007) and cleaned reactions using a 125 mM ethylenediaminetetraacetic acid–ethanol protocol. We conducted automated sequencing of complimentary strands using the Applied Biosystems 3110 DNA sequencer of the Molecular Biological Facility at the University of New Mexico.
We edited and aligned sequences in sequencher 4.8 (Genecodes, Ann Arbor, Michigan) and verified visually. We translated protein coding sequences for Cytb and ApoB genes to amino acids and examined them for internal stop codons, rates of transition (Ti)–transversion (Tv) changes, and relative 1st, 2nd, and 3rd position changes in codons that might indicate a pseudogene. We compared complementary strands of DNA and deposited contiguous sequences in GenBank (Appendix I). We inferred alleles of nuclear heterozygotes using phase (Stephens et al. 2001), which implements a Bayesian method for reconstructing haplotypes from nuclear sequences that include multiple heterozygous base sites within individuals. phase was run to accept results with a probability > 90%, using a burn-in of 10,000 iterations and a run length of 10,000 iterations (McCormack et al. 2011). Among 5 runs, we retained results from the run with best goodness-of-fit to an approximate coalescent model, resulting in 2 nuclear haplotype sequences or alleles per individual. We used only 1 allele per individual (chosen randomly) for analyses.
Gene tree relationships
We estimated independent genealogies for all loci. To retain maximum genetic diversity in the Cytb data set, we included incomplete sequences and did not assign haplotypes. All nuclear genealogies also utilized all available sequences. For the Cytb genealogy (n = 168), we removed 37 previously published sequences to allow for ease of visual interpretation. All excluded sequences represented duplicate haplotypes that we included in subsequent population demographic analyses. Although we acknowledge the potential for bias using all individuals, estimation of species tree relationships benefit from inclusion of all sequences (McCormack et al. 2011). We determined DNA substitution models for each locus using MrModeltest version 2.3 (Nylander 2004) under the Akaike information criterion. We performed Markov chain Monte Carlo searches in MrBayes version 3.1 (Ronquist and Huelsenbeck 2003). Each of 3 independent runs computed 10 million generations, sampling every 1,000 generations, with 5 Markov chains and default heating values, the first 1,000 trees discarded as burn-in, and the short-branch method (Marshall 2010). Stationarity of Markov chain Monte Carlo runs was assessed in Tracer version 1.4 (Rambaut and Drummond 2007). Phylograms were midpoint rooted due to uncertain systematic relationships within the vagrans complex (no unambiguous sister taxon to the North American water shrews), and were visualized with posterior probabilities in FigTree version 1.2.2 (Rambaut 2009).
Clade demographics
Genetic diversity analyses for the Cytb locus used all available North American water shrew sequences, with diversity and demographic analyses performed separately for major group assignments. We calculated average pairwise and corrected sequence divergence among groups for 871 bp of Cytb (positions 130–1,001), the largest length available for all groups. In addition, we tested hypotheses of expansion versus stability across latitude in the Cordilleran group. We grouped samples from areas previously covered by continental ice sheets (Cordilleran North) separate from samples from southern unglaciated areas (Cordilleran South), except for southwestern localities representing distinct genetic lineages. We calculated summary statistics in DnaSP (Librado and Rozas 2009) and assessed genetic diversity including number of haplotypes (h), haplotype diversity (Hd), nucleotide diversity (π), and pairwise sequence divergence. To estimate relative population sizes, we calculated theta (θ), using LAMARC version 2.1.6 (Kuhner 2006). For tests of demographic expansion, we used DnaSP to calculate Tajima's D (Tajima 1989), Fu's FS (Fu 1997), and R2 (Ramos-Onsins and Rozas 2002) and assessed significance with 10,000 coalescent simulations. We also implemented the Markov chain Monte Carlo approach using the Metropolis–Hasting algorithm (Kuhner et al. 1998) in LAMARC to investigate the population growth parameter (g) for each locus and clade. For each run within a maximum-likelihood framework and a starting value of g = 1, we used 1,000 short chains (sampling increments of 10; 500 trees sampled per chain; 1,000 trees discarded), 10 long chains (sampling increments of 10; 20,000 trees sampled per chain; 1,000 trees discarded), random starting seed, and 2 simultaneous searches with heating temperature of 1.2. We used the F84 substitution model with empirical base frequencies and a Ti/Tv ratio = 15. Results of 2 independent runs with different random starting seeds for each test were assessed for convergence in tracer (effective sample size ≥ 300) and to confirm that different runs produced equivalent results. To avoid potential upward bias, we calculated the standard deviation from each mean value of g and inferred significant population growth if g > 3 SD (Lessa et al. 2003).
To calculate population size change through time, we used BEAST version 1.6.1 (Drummond and Rambaut 2007) using the coalescent Bayesian Skyline Plot tree prior. Results of Bayes factor tests to determine the best clock model indicated no significant difference between strict and relaxed clocks. Subsequently we utilized a strict clock considering no analyses were found to deviate from a constant evolutionary rate. Parameters were set in BEAUti, part of the BEAST software package. We determined the appropriate model of nucleotide substitution for each lineage using MrModeltest version 2.3 under the Akaike information criterion. We ran Markov chain Monte Carlo analyses (length of chain = 100 M logging trees every 1,000) with a mutation rate of 5.5% per million years, previously estimated for another species of Sorex (Hope et al. 2010). We assume that closely related taxa share similar substitution rates through time (Kumar and Subramanian 2002).
Species tree relationships and divergence dating
Multiple independent genes are critical for inferring species relationships (Maddison 1997; Edwards 2009) because gene duplication, introgression, sorting of ancestral polymorphism, or rapid fixation of linked loci through selective sweeps can be problematic with single genealogies (Carstens and Knowles 2007). We estimated species tree relationships with *BEAST (Heled and Drummond 2010) using a method that involved coalescent coestimation of genealogies from all 4 loci embedded within a corresponding species tree topology and implemented in the program BEAST. This method uses a Bayesian Markov chain Monte Carlo multispecies coalescent technique. Because *BEAST accommodates different numbers of gene copies for each taxon (Heled and Drummond 2010), we included all samples. We unlinked data sets for each locus across all partitions, assigned substitution models by locus, and applied empirical base frequencies. We designated samples to major groups of North American water shrews or to the various outgroup taxa, excluding 3 individuals with incongruent matrilineal and nuclear signatures between Cordilleran and Boreal groups from analysis. We applied a Cytb mutation rate of 5.5% per million years to that genealogy only, leaving rates based on other loci to be estimated. From Bayes factor tests of preliminary runs, we applied a relaxed clock (uncorrelated lognormal) to Cytb as it consistently produced lowest −lnL scores and because evolutionary rates vary at different depths within a phylogeny (Ho et al. 2011). We assigned a strict clock to nuclear loci to reflect lower divergence over the evolutionary time frame investigated and chose a Yule tree prior with piecewise linear and constant root. We assigned proper ploidy to each locus. We found that default values for priors and operators produced acceptable effective sample size (ESS) values from preliminary runs. Two independent runs used Markov chain Monte Carlo chains of 100 million generations sampling trees every 10,000. We annotated species tree files in TreeAnnotator (BEAST package) and visualized the species tree in FigTree. We retrieved all node values, ages (heights), and associated 95% confidence intervals from tree files.
To investigate differences based on all loci versus only the mitochondrial locus, we ran further *BEAST analyses using only the Cytb data set or only the 3 nuclear loci. Although rapid sorting of mitochondrial genealogies may detect recent diversification, more slowly evolving nuclear loci may control for potential homoplasy deeper in the phylogeny. In addition, estimates from mitochondrial genes alone may better account for the actual timing of gene divergence that predates species divergence (Heled and Drummond 2010). We repeated species tree estimation twice using all loci but with different levels of resolution for lineages (assignments based on the Cytb phylogeny) to test for consistency of divergence estimates using different lineage constraints. Lower resolution considered main groups of North American water shrews as originally described (Coastal, V.I., Cordilleran, Boreal, and Eastern), higher resolution analysis considered additional clade assignments by splitting the Cordilleran group into 3 subgroups (Main, Southwest, and White Mtns.), V.I. into 2 lineages (1 and 2 reflecting possibly multiple colonizations of Vancouver Island), and Eastern into 2 subgroups representing the subspecies albibarbis and gloveralleni.
Results
Gene tree relationships
Independent genealogies for the 4 separate loci ranged from well sorted to unresolved among the North American water shrew lineages (Fig. 2; Supporting Information S1–S3, DOI: 10.1644/13-MAMM-A-196.S1, 10.1644/13-MAMM-A-196.S2, and 10.1644/13-MAMM-A-196.S3). Both Cytb and ApoB loci indicated a lack of monophyly for North American water shrews that were paraphyletic with respect to other species of the vagrans complex. Both also indicated that Cordilleran and V.I. groups form a single intermixed lineage. Southwest and White Mtns. lineages of the Cordilleran group are each well supported based on the Cytb locus and also exhibited significant divergence within the ApoB genealogy. The Cordilleran group also exhibited regional geographic structure consistent with a previous mtDNA assessment (Himes and Kenagy 2010) indicating distinct Cascades, Northern Rockies, and Olympic Peninsula subclades, as well as 2 well-supported subclades from V.I. Coastal and Cordilleran groups represent reciprocally monophyletic sister lineages, as do Boreal and Eastern groups (Cytb). Only IFN1 indicated monophyly of all North American water shrews. MYH2 largely is unresolved within North American water shrews, although outgroup taxa exhibit greater divergence. Specimens belonging to different subspecies within the Eastern group clearly are not resolved over any locus. Three individuals (RBCM:020008, RBCM:020016, and RBCM:020017 [Fig. 1—localities 44, 46, and 45, respectively]) exhibited incongruent matrilineal (Boreal [Fig. 2]) and nuclear (Cordilleran [Supporting Information S1]) assignments. A well-supported clade including S. vagrans and S. ornatus consistently exhibited reciprocal monophyly with regard to other species of the vagrans complex over all loci except MYH2.
Fig. 2.
Bayesian phylogeny of North American water shrews and related species representing a midpoint-rooted genealogy based on mitochondrial DNA cytochrome-b (Cytb) sequences (639–1,140 bp). Posterior probabilities for nodes supported ≥ 0.95 are indicated by black dots. Probability values < 0.90 are omitted. Asterisk (*) indicates incongruent matrilineal individuals compared with the nuclear apolipoprotein B exon (ApoB) genealogy (RBCM:020008, RBCM:020016, and RBCM:020017; localities 44, 46, and 45, respectively). Locality numbers correspond to Appendix I. V.I. clades 1 and 2 refer to samples used for species-tree analysis.
Genetic diversity and clade demographics
The Eastern group exhibited highest Cytb nucleotide and haplotype diversity, followed by the Cordilleran group (Table 1). Lowest diversity was in Boreal followed by V.I. groups, whereas the Coastal group exhibited intermediate diversity. Relative values of theta for different groups, reflecting contemporary effective population size, were comparable to nucleotide diversity. Highest effective population sizes were in Eastern and Cordilleran groups and lowest were in V.I. and Boreal groups. Expansion statistics were significant for Cordilleran and Eastern groups and nonsignificant for V.I. and Boreal groups. The Coastal group only exhibited significant population expansion based on the g-statistic.
Table 1.
Genetic diversity indexes and growth estimates for all groups of North American water shrews based on mitochondrial DNA cytochrome-b sequence data. n = sample size; L = length of sequence; h = number of haplotypes; Hd = haplotype diversity; π = nucleotide diversity; θ = theta reflecting effective population size; D = Tajima's D (with associated P-value); FS = Fu's FS (with associated P-value); R2 = Ramos-Onsin and Rozas's R2 (with associated P-value); g = growth statistic (with associated standard deviation); SD(g) > 0 = number of standard deviations above zero of the growth statistic. Values in boldface type represent significant growth statistics at the P ≤ 0.05 level, or for ≥ 3 SD(g) > 0.
Average pairwise sequence divergence (based on Cytb) between taxa within the vagrans complex ranged from 0.56% to 9.09% (Table 2). Cordilleran and V.I. groups were < 1% divergent. Coastal shrews were ∼3% divergent from Cordilleran and V.I. water shrews and ∼6% divergent from Boreal and Eastern water shrews. These latter 2 groups were ∼3% divergent from each other, but were each ∼5–7% divergent from other North American water shrews. Pairwise divergence between S. monticolus, S. pacificus, and S. sonomae and divergence of each from all North American water shrews were comparable to levels of divergence between North American water shrew groups. S. ornatus and S. vagrans exhibited highest divergence values except from each other.
Table 2.
Percent sequence divergence (upper right of diagonal) and corrected sequence divergence (lower left of diagonal) between groups of North American water shrews. Coastal group represents Sorex bendirii and other groups (2–9) represent lineages within S. palustris.
Considering an estimated Cytb mutation rate of 5.5% per million years, timing of demographic population size changes varied among North American water shrew lineages (Fig. 3; Supporting Information S4, DOI: 10.1644/13-MAMM-A-196.S4). The Coastal group did not exhibit significant size change through time considering broad confidence limits. Extreme population growth occurred in Cordilleran (North), Eastern, and Cordilleran (South) groups, with most growth occurring before the Last Glacial Maximum, coincident with interstadials of the mid- to late Wisconsinan glacial period (75–25 kya), then levelling to the present. Slight (nonsignificant) population growth of both V.I. and Boreal groups was dated to post–Last Glacial Maximum.
Fig. 3.
Bayesian skyline plots for lineages of North American water shrews created in BEAST from the mitochondrial DNA cytochrome-b (Cytb) data set. Lines indicate change in effective population size through time. Asterisks (*) indicate that population size change does not breach 95% confidence interval (see Supporting Information S4 for individual group Bayesian skyline plots). Bayesian skyline plots extend right to left from past to present scaled in thousands of years. Vertical axis represents τ as a function of generation time and effective population size. Vertical dashed line within the plot indicates the Last Glacial Maximum placed at 21 kya.
Species tree relationships and divergence dating
Phylogenetic relationships based on all 4 loci supported major North American water shrew lineages as recognized from previous studies (O'Neill et al. 2005; Himes and Kenagy 2010; Mycroft et al. 2011). But there was no posterior support for monophyly of all North American water shrews despite the consensus topology (Fig. 4) except for moderately high support when considering only the 3 nuclear loci (Supporting Information S5, DOI: 10.1644/13-MAMM-A-196.S5). Consistent monophyly of an Eastern–Boreal clade was well supported, as was monophyly of a Cordilleran–V.I. clade. When Southwest and White Mtns. clades of the Cordilleran group were constrained as independent lineages, they were rendered basal to the combined Main–V.I. clade (Fig. 4b). With an understanding that coalescence times generally predate lineage splitting and should be considered as maximum divergence estimates, coalescence of the Cordilleran–V.I. clade was consistently estimated at 0.05 mya, although Southwest and White Mtns. lineages were dated to 0.11 mya and 0.21 mya, respectively. The Eastern–Boreal clade coalesces at 0.29 mya. Estimates for initial divergence within the vagrans complex are ∼1 mya.
Fig. 4.
Species-tree estimation of North American water shrews and related species for a) main clade, and b) additional clade assignments based on 4 independent loci including mitochondrial DNA cytochrome-b (Cytb) gene, and regions from nuclear apolipoprotein B, interferon 1, and myosin heavy-chain 2 genes. Phylogeny estimation was performed in *BEAST, providing both posterior probability nodal support values (bold, left of node) and coalescence times (millions of years; italics, right of node). Shaded bars around nodes indicate 95% confidence intervals for coalescence estimates. Topologies are presented as ultrametric and proportional to an evolutionary timeline (bottom) that extends from the present (right) to the past (left). Asterisks (*) highlight low posterior support for basal nodes of the vagrans complex, indicating a lack of monophyly for North American water shrews.
Discussion
Widespread circumboreal forests constitute the largest biome in the world and support a rich mammalian fauna. At a regional scale, constituent species often share aspects of their evolutionary histories (Arbogast and Kenagy 2001) because Quaternary environments predictably influenced diversification processes (Hewitt 1996). As more species are examined, generalized patterns are emerging, with differences among them providing insight into how species-specific attributes, such as vagility, habitat affinities, competition, and life-history differences may effect species' responses to environmental change. More vagile species that are widely distributed and easily maintain population connectivity across a landscape may contain relatively few distinct lineages (Dawson et al. 2014), whereas species less able to cross persistent or even intermittent barriers may exhibit finer-scale phylogeographic structure (Demboski and Cook 2001; O'Neill et al. 2005; Hafner and Smith 2010; Himes and Kenagy 2010; Hope et al. 2011, 2012). Our investigation of North American water shrews using multiple independent loci revealed a phylogeographic history that has culminated in multiple distinct species, making these semiaquatic forest dwellers a valuable comparative model for other codistributed taxa.
Concerning systematic relationships, independent genealogies (except for IFN1) are inconsistent with a single origin of North American water shrews with respect to other species of the vagrans complex and there is minimal nodal support for a North American water shrew clade from the species tree analyses. None of the insular North Pacific Coast populations, including S. p. brooksi from Vancouver Island and 2 island samples of S. p. navigator from southeastern Alaska, is significantly differentiated from the mainland samples of the Cordilleran group, indicating no signature of divergence due to coastal Quaternary refugia. However, the Cordilleran (combined with V.I.), Coastal, Boreal, and Eastern groups all exhibit moderately high genetic divergence from each other (Table 2). Support for 4 species of North American water shrews therefore is reflected by polyphyly of all North American water shrews with respect to other species of the vagrans complex considering both mitochondrial and nuclear loci; reciprocal monophyly of all lineages representing nominal species, each with branch lengths from species tree reconstructions that reflect prolonged divergence spanning several late-Quaternary glacial cycles (support for the phylogenetic species concept—Cracraft 1983); substantial geographic structure although with potential zones of contact between lineages; and apparent genetic isolation of lineages despite contact (biological species concept), except for 3 specimens with mitochondrial incongruence that may reflect lingering unsorted ancestral haplotypes, ancient admixture, or possibly a lack of complete reproductive isolation between the 2 most divergent North American water shrew lineages.
We therefore recommend recognition of 4 species of North American water shrews represented by Coastal (Sorex bendirii [Merriam] 1884; marsh shrew), Cordilleran and V.I. combined (Sorex navigator [Baird] 1858; western water shrew), Boreal (Sorex palustris Richardson, 1828; American water shrew), and Eastern (Sorex albibarbis [Cope] 1862; eastern water shrew) groups. Finally, we highlight the tentative specific status of Sorex alaskanus Merriam, 1900, pending genetic analysis. As reflected by mitochondrial divergence from previous studies, the current sampling also indicated levels of mtDNA divergence in the range seen between other recognized species of Sorex (Hope et al. 2012).
Our data support O'Neill et al. (2005), who recognized S. navigator as distinct from North American water shrews occurring east of the western Cordillera (Boreal and Eastern). Jackson (1928) and Soper (1964) noted that specimens of the Cordilleran group differed from the Boreal group in dorsal pelage color, cranial morphology, and dentition. Despite limited incongruent mitochondrial haplotypes, genetic distance between these taxa is the greatest among all North American water shrew comparisons (Table 2; Mycroft et al. 2011), indicating accumulated divergence over ∼0.5 × 106 years (Fig. 4). Possible parapatry through Alberta and northeastern British Columbia between Cordilleran and Boreal groups may represent contact following extended isolation. Alternatively, repeated contact with limited introgression may have occurred over multiple glacial cycles (i.e., reinforcing selection). Either scenario lends support for specific status between S. navigator and S. palustris. Further sampling including voucher specimens of Cordilleran and Boreal groups would help to characterize the extent and dynamics of contact. We also maintain synonymy of the subspecies S. p. hydrobadistes within S. palustris considering geographic distribution, although it occurs in a poorly sampled region between both Boreal and Eastern groups and further genetic analyses should confirm this designation.
Within the Cordilleran group (S. navigator), we report 2 distinct genetic lineages in mountains of the southwestern United States, with water shrews from the White Mtns. of Arizona the most divergent from other Cordilleran water shrews. We retain these lineages as S. navigator considering our limited sampling and no reciprocal monophyly with the remainder of the Cordilleran group. However, the restricted and isolated geographic ranges of these lineages in southwestern montane habitats indicate the need for further investigation of distributional limits and population status given few available samples and because of current warming and drought-prone climate trends in this region (Smith 1993; MacDonald 2010). Minimally, these lineages represent distinct evolutionarily significant units (Moritz 1994) and may warrant subspecific status. Although genetically similar to, and polyphyletic with respect to the remainder of the Cordilleran group, we retain the V.I. group as a subspecies, classifying it as S. navigator brooksi, pending expanded morphological and genetic analyses.
Our results and the analysis by Mycroft et al. (2011) support recognition of the Eastern group as a distinct species S. albibarbis (including subspecies S. p. albibarbis, S. p. gloveralleni, and S. p. labradorensis) because of reciprocal monophyly of Eastern and Boreal groups, no evidence of contemporary gene flow from the small sample included, and moderate genetic divergence (3.22%). Based on geographic distribution of subspecies S. p. punctulatus and S. p. turneri we have synonimized them with S. albibarbis, although further genetic sampling should confirm their specific status. Although genetic divergence of Cytb is modest compared with other recognized species-level divergences, it is higher between Eastern and Boreal groups than between Coastal and Cordilleran groups of water shrews, the latter 2 representing both genetically isolated and morphologically distinct species.
Our data support the recognition by O'Neill et al. (2005) of S. bendirii as distinct from other North American water shrews due to genetic differentiation and lack of evidence of gene flow despite parapatry (Fig. 1). In addition, S. bendirii is morphologically (dark venter and indistinct toe fringes versus silver-gray venter and distinct toe fringes in S. palustris) and ecologically (broadly associated with wet habitats versus a more strict association with running streams for S. palustris) distinctive, although some overlap exists (Pattie 1973; van Zyll de Jong 1983; Beneski and Stinson 1987). Merriam (1884) and Jackson (1928) observed morphological convergence in some characters but concluded there was no evidence for morphological intergradation.
Clarification of evolutionary relationships within the vagrans complex has been problematic because of morphological similarity among several species (Findley 1955; Carraway 1990) and apparent disagreement between molecular and morphological evidence (Demboski and Cook 2001; Shafer and Stewart 2007). Carraway (1990) noted that the “Sorex vagrans complex” is a misnomer because S. vagrans is not the most ancestral member of the group. Analyses of multiple loci indicate a close sister association of S. vagrans and S. ornatus with all North American water shrews and other recognized members of the vagrans complex included here, but we retain reference to the vagrans complex because of the lack of an unambiguous ancestral lineage (Fig. 4).
Although species tree estimation indicates a well-supported early dichotomy within the vagrans complex dated to ∼1 mya, there is consistent lack of nodal support between North American water shrews and other members of the complex. It is possible that differentiation within the vagrans complex resulted from initial fragmentation (hard polytomy) around 0.5–0.7 mya with or without ongoing gene flow through time considering their close geographic proximity. Ancient hybridization between S. palustris and S. monticolus has been postulated as accounting for lack of monophyly of North American water shrews (Demboski and Cook 2001; Shafer and Stewart 2007) and may indicate a speciation-with-gene-flow mechanism for diversification within this shrew complex (Feder et al. 2012). Frequency of divergence events dated through coalescent species trees are generally coincident with periodicity of glacial cycling (∼100,000 years) since the middle Quaternary, indicating that climate cycling strongly influenced speciation within this diverse group of shrews, possibly due to repeated episodes of allopatry and divergence as seen for other taxa (Galbreath et al. 2010; Hope et al. 2012; Dawson et al. 2014) and possibly with intermittent hybridization. Future genomic investigation would provide a more detailed assessment of divergence and levels of gene flow. At some time earlier than 0.3 mya, ancestors of Boreal and Eastern groups likely were isolated farther east in North America and diverged in allopatry from all other species of the vagrans complex. Twice (at 0.21 mya, pre-Illinoian interglacial; and at 0.11 mya, Sangamon interglacial) Cordilleran water shrews were isolated in the southwestern United States, diverging subsequently. Most recently water shrews of the Cordilleran group colonized islands of the North Pacific Coast including both Vancouver Island and a few islands of the Alexander Archipelago of southeastern Alaska as the Cordilleran continental ice sheet receded following the Wisconsinan glacial period (< 50 kya). The presence of 2 distinct mtDNA subclades within the V.I. group and other samples negligibly divergent from mainland specimens of the Cordilleran group suggests that Vancouver Island was colonized on multiple occasions post–Last Glacial Maximum, or alternatively that both unique and ancestral haplotypes are retained within the V.I. group following a single colonization. Most parsimoniously, colonization of Vancouver Island would have been via a stepping-stone route through the San Juan Islands, a recognized filter bridge (Wilson et al. 2009) or possibly from the Olympic Peninsula to the south when straits were narrower following the Last Glacial Maximum. Insular populations in southeastern Alaska show no genetic divergence from mainland Cordilleran specimens. Only a few large coastal islands in North America support populations of North American water shrews, suggesting relatively low dispersal ability across marine barriers or relatively large habitat requirements to sustain viable populations, or both.
Fossil evidence is relatively abundant for North American water shrews (Fig. 1; Kurtén and Anderson 1980). Himes and Kenagy (2010) concluded that the Cordilleran group (not including samples from the southwestern United States) reached their current distribution by dramatic range expansion from a single Last Glacial refugium south of the Cordilleran ice sheet. Our results support this conclusion although fossils in the southern Rockies may represent other isolated refugia that correspond to distinct southwestern United States clades. Support for regional Cascade and Northern Rockies clades within the Cordilleran group indicates a history of regional isolation and reconnection in response to recent climate cycles that is paralleled in other montane species through the western United States (e.g., Galbreath et al. 2010; Hafner and Smith 2010). In addition, although fossils from mesic forested regions such as the North Pacific Coast are rare, a distinct genetic lineage centered on the Olympic Peninsula (Figs. 1 and 2; Himes and Kenagy 2010) conceivably represents a refugial population from this once periglacial area. Other Wisconsinan fossils from Kansas eastward to Pennsylvania may reflect refugia inhabited by Boreal and Eastern groups during this glacial phase. Numerous Wisconsinan fossil localities have yielded evidence of multiple forest-associated species much farther south than the present distribution of the boreal biome (Kurtén and Anderson 1980). Inclusion of ancient DNA in future comparative phylogeographic investigation may further clarify location of subrefugia for different species and regional phylogroups of North American water shrews.
Although dates of major nodes within the species tree appear robust to different levels of lineage resolution based on group designations, we caution interpretation of coalescence times from species trees based on only a few differentially informative loci (McCormack and Faircloth 2013). For instance, we dated coalescence of S. cinereus and S. haydeni to 0.81 mya although the same node from species tree estimation of the cinereus complex of shrews by Hope et al. (2012—see Fig. 3) was dated earlier. However, coalescence of S. palustris and S. monticolus from Hope et al. (2012—see Fig. 3) was dated to 0.55 mya, comparable to dates from current species trees (Fig. 4). We suspect that dating discrepancies between the 2 analyses may reflect artifacts of sampling relatively few loci, and few individuals for outgroup taxa, as well as different evolutionary histories between the vagrans and cinereus species complexes.
Regarding sampling, different outgroups and loci have been included in phylogeographic analyses of shrews (Hope et al. 2012). Additional loci can affect estimated demographic and phylogenetic parameters, depending on information content, mutation rate, and particularly the total number of loci and individuals (Maddison and Knowles 2006; Carstens and Knowles 2007). Using too few loci may alter coalescence estimates (Heled and Drummond 2010), although too many loci also may detrimentally influence coalescent inference (McCormack and Faircloth 2013). Only 1 S. cinereus and 1 S. haydeni were included in the species trees presented herein, perhaps inflating the estimation of fixed differences between the 2 samples over all loci. In Hope et al. (2012), multiple samples represented each species within the cinereus complex and consequently, incomplete lineage sorting or contemporary gene flow between lineages, or both, resulted in much shallower coalescence (Carstens and Knowles 2007).
Finally, we illustrated that lineage designations significantly influenced coalescence estimates. Species tree estimation for North American water shrews, grouping all Cordilleran lineages (including the 2 southwestern United States lineages) and all V.I. specimens together, resulted in a coalescent estimate for the Cordilleran–V.I. clade at 0.05 mya. However, when lineage designations were included to recognize additional genetic diversity of the southwestern United States populations, coalescence of the Cordilleran–V.I. clade remained at 0.05 mya but coalescences of additional Southwest and White Mtns. lineages were much deeper. Subsuming distinct lineages represented by very few individuals within a much larger clade therefore may obscure signals of earlier divergence events.
A general lack of population size change for the Coastal group is consistent with persistence in western coastal forests over extended periods. This area was never covered by ice sheets and mesic forests prevailed during both glacial and interglacial periods. The V.I. group also exhibits population stability, perhaps reflecting confinement within a small area and potentially multiple colonizations from nearby mainland coastal forests. Although there is no signal of a refugial lineage along the North Pacific Coast that persisted through the last glacial, Vancouver Island may have been colonized soon after the Last Glacial Maximum when coastal areas were exposed from under continental ice (Nagorsen and Keddie 2000; Brown and Hebda 2003). The timing of significant expansion of both “North” and “South” Cordilleran groups indicates not only recolonization of previously ice-covered areas, but also regional expansion south of continental ice sheets (Himes and Kenagy 2010).
Following glacial retreat, Cordilleran, Boreal, and Eastern North American water shrews expanded northward but there is evidence that this was not a rapid pioneer advance to follow northward colonization of forests. First, significant population growth as recorded in the genetic legacies of lineages that now occupy high latitudes occurred prior to ice-sheet recession. In addition, the Boreal group exhibits no evidence of significant population expansion or growth; this may instead reflect a signature of slower phalanx expansion in this region (Hewitt 1996; Yan et al. 2012). Modern range limits of North American water shrews do not reach the northern extent of available forested habitats at high latitude (Banfield 1977; van Zyll de Jong 1983; Cook et al. 1997; MacDonald and Cook 2009). This delayed northward movement contrasts with numerous other forest-associated species that rapidly expanded to beyond the northern extent of boreal biomes in North America following the Last Glacial Maximum (Demboski and Cook 2001; Jaramillo-Correa et al. 2004; Hope et al. 2012). We suggest that the semiaquatic habit of North American water shrews may entail a northern limit for these small mammals that is coincident with streams and lakes that remain at least partially unfrozen through winter months. Forests may well provide an additional overstory buffer for stream access (and hence aquatic invertebrate food source) through cold periods (Beneski and Stinson 1987). Such an association has strong implications for continued northward expansion of North American water shrews in response to a warming Arctic.
There are a number of management implications considering distinct lineages. Our multilocus approach for testing species limits among major North American water shrew lineages confirmed the species status of S. bendirii. British Columbian populations of this species are listed as endangered under Canada's Species At Risk Act because of habitat loss from urban development (COSEWIC 2006). However, additional surveys and voucher collection north and east of the known Canadian range of S. bendirii are required to determine its geographic limits and degree of parapatry with S. navigator. Other potentially distinct taxa (Hall 1981) should also be investigated. We have shown that both the V.I. subspecies (S. n. brooksi) and samples collected from the Alexander Archipelago are only minimally divergent from mainland S. navigator. Evidence for a Last Glacial Maximum refugium on Vancouver Island is weak (Nagorsen and Cardini 2009) and though coalescence estimates for V.I. clades marginally predate the Last Glacial Maximum (Fig. 4), minimally divergent V.I. lineages indicate more recent isolation. In addition, insularity may stimulate rapid morphological differentiation (Cardini et al. 2009) and elevated mutation rates (Smith and Klicka 2013) and despite genetic similarity between North Pacific Coast and Cordilleran specimens, S. n. brooksi and Alexander Archipelago populations represent insular isolates on unique evolutionary trajectories. These populations persist on islands subject to increasing anthropogenic as well as stochastic processes, therefore taxon status as representative of unique regional faunas may warrant consideration in future management initiatives (Craig 2004). Southwest and White Mtns. lineages of S. navigator are montane isolates within sky-island forested habitat where increasing fire hazards and forest diseases are a concern for persistence of endemic relictual lineages (MacDonald 2010). In addition, close associations of North American water shrews with increasingly fragile and isolated riparian zones in the southwestern United States may compound risk of extirpation (Smith 1993).
More comprehensive genetic sampling, particularly from central and eastern North America, representing all subspecies will facilitate more rigorous identification of the geographic limits of each North American water shrew taxon (and lineage), and characterization of the location and extent of contact zones (Fig. 1). Ongoing morphological assessments will help to corroborate or refute current systematic relationships now based on genetic data, and to develop a set of diagnostic morphological traits for each species. Further, a comparative phylogeographic perspective for North American water shrews of the vagrans complex would present an excellent opportunity to elucidate fundamental speciation mechanisms. Rapid evolution within the cinereus complex of shrews is a consequence of repeated allopatry followed by expansion of lineages at different times, in accordance with different ecological associations (Hope et al. 2012). Although comparably diverse, the vagrans complex instead exhibits extensive sympatry between most species, including S. bendirii and S. navigator. Only the eastern species S. palustris and S. albibarbis are allopatric with respect to the remainder of the vagrans complex. With the growing efficacy of genomic methods (e.g., Feder et al. 2012), testing evolutionary hypotheses related to timing and sequence of diversification should be productive, particularly considering potential for ongoing gene flow between the 2 most genetically divergent North American water shrew lineages (Cordilleran and Boreal).
This study highlights the utility of wide-ranging species complexes for evolutionary studies. Strong habitat associations and distinct phylogeographic structure within North American water shrews constitutes an excellent system for investigating macroecological and macroevolutionary processes at the interface of terrestrial and aquatic environments.
Acknowledgments
We thank T. Jung, J. Whitman, S. O. MacDonald, C. Ramotnik, J. Malaney, and J. Lackey for assistance with sample collection. Tissues were supplied from the Museum of Southwestern Biology, Royal British Columbia Museum, and the University of Alaska Museum of the North. Funding was provided by the National Science Foundation (NSF0415668, 0731350, and 1258010) and the United States Geological Survey's Alaska Regional Director's DOI on the Landscape initiative. Statistical analyses were facilitated by the University of Alaska, Fairbanks, Life Science Informatics Portal, http://biotech.inbre.alaska.edu. University of Alaska, Fairbanks, Life Science Informatics as a core research resource is supported by grant RR016466 from the National Center for Research Resources, a component of the National Institutes of Health. Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the United States Government. We thank D. J. Hafner for insightful discussion and comments on an earlier version, and J. Pearce and 2 anonymous peers for reviewing the manuscript.
Supporting Information
Supporting Information S1.—Bayesian phylogeny of North American water shrews and related species representing a midpoint-rooted genealogy based on nuclear apolipoprotein B exon sequences (434 bp). Posterior probabilities for major nodes are illustrated. Asterisks (*) indicate individuals with incongruent matrilineal and nuclear lineage associations. Locality numbers correspond to Appendix I.
Found at DOI: 10.1644/13-MAMM-A-196.S1 (10.1644_13-MAMM-A-196.S1.eps)
Supporting Information S2.—Bayesian phylogeny of North American water shrews and related species representing a midpoint-rooted genealogy based on nuclear interferon 1 intron sequences (282 bp). Posterior probabilities for major nodes are illustrated. Asterisks (*) indicate individuals with incongruent matrilineal and nuclear lineage associations. Locality numbers correspond to Appendix I.
Found at DOI: 10.1644/13-MAMM-A-196.S2 (10.1644_13-MAMM-A-196.S2.eps)
Supporting Information S3.—Bayesian phylogeny of North American water shrews and related species representing a midpoint-rooted genealogy based on nuclear myosin heavy-chain 2 intron sequences (286 bp). Posterior probabilities for major nodes are illustrated. Asterisks (*) indicate individuals with incongruent matrilineal and nuclear lineage associations. Locality numbers correspond to Appendix I.
Found at DOI: 10.1644/13-MAMM-A-196.S3 (10.1644_13-MAMM-A-196.S3.eps)
Supporting Information S4.—Bayesian skyline plots for lineages of North American water shrews created in BEAST from the mitochondrial DNA cytochrome-b (Cytb) data set. Middle lines indicate change in effective population size through time accompanied by 95% confidence intervals (shaded interval). Bayesian skyline plots extend right to left from past to present scaled in thousands of years. Vertical axis represents τ as a function of generation time and effective population size.
Found at DOI: 10.1644/13-MAMM-A-196.S4 (10.1644_13-MAMM-A-196.S4.eps)
Supporting Information S5.—Species tree estimation of North American water shrews and related species based on a) 3 independent nuclear loci including regions from nuclear apolipoprotein B, interferon 1, and myosin heavy-chain 2 genes, and mitochondrial DNA cytochrome-b (Cytb) gene. Phylogeny estimation was performed in *BEAST, providing both posterior probability nodal support values (bold, left of node) and coalescence times (millions of years; italics, right of node). Shaded bars around nodes indicate 95% confidence intervals for coalescence estimates. Topologies are presented as ultrametric and proportional to an evolutionary timeline (bottom) that extends from the present (right) to the past (left).
Found at DOI: 10.1644/13-MAMM-A-196.S5 (10.1644_13-MAMM-A-196.S5.eps)
Literature Cited
Appendices
Appendix I
Specimens examined are listed by scientific name (reflecting updated taxonomy), specific locality number (in boldface type, corresponding to Fig. 1), specific locality, specimen number (museum catalog number), and corresponding GenBank accession numbers (in order: mitochondrial DNA cytochrome-b [Cytb], nuclear DNA apolipoprotein B, interferon 1, and myosin heavy-chain 2; single accessions for specimens reflect only Cytb sequences). Acronyms for museum archives are DMNS = Denver Museum of Nature & Science; MSB = Museum of Southwestern Biology; RBCM = Royal British Columbia Museum; ROM = Royal Ontario Museum; UAM = University of Alaska Museum of the North, Fairbanks; UWBM = University of Washington, Burke Museum of Natural History and Culture. Sequence accessions previously published and retrieved from GenBank for analysis in the current study are also listed by original citation. NA = not applicable.
Sorex navigator (brooksi).—(1) Canada, British Columbia: Vancouver Island, Kay Creek, 50.1486°, −125.7639°: RBCM019738 (KF302816); Miller Creek, 50.0025°, −125.5667°: RBCM019865 (KF302818, KF302890, KF302969, KF303047); North Memekay River, 50.1856°, −125.9408°: RBCM019983 (KF302821); (2) Canada, British Columbia: Vancouver Island, Woodhus Creek, 49.8853°, −125.2958°: RBCM019721 (KF302807); Black Creek, 49.8703°, −125.2569°: RBCM019722 (KF302808); RBCM019723 (KF302809, KF302882, KF302961, KF303039); 49.8608°, −125.2203°: RBCM019867 (KF302819); 49.8592°, −125.2439°: RBCM020309 (KF302824); Tsolum Creek, 49.8161°, −125.2136°: RBCM019724 (NA, KF302883, KF302962, KF303040); (3) Canada, British Columbia: Vancouver Island, Lowry Lake, 49.4000°, −125.1333°: RBCM019496 (AY954928); RBCM019497 (AY954929); RBCM019498 (AY954930); Dove Creek, 49.7244°, −125.0978°: RBCM020254 (KF302822); Morrison Creek, 49.6500°, −125.0294°: RBCM020256 (KF302823); (4) Canada, British Columbia: Vancouver Island, Hamilton Creek, 49.3167°, −124.4333°: RBCM019493 (AY954926); RBCM019495 (AY954927); Black Brook, 49.3956°, −124.6344°: RBCM019730 (KF302810, KF302884, KF302963, KF303041); Morrison Creek, 49.2692°, −124.3544°: RBCM019731 (KF302811, KF302885, KF302964, KF303042); (5) Canada, British Columbia: Vancouver Island, Miller Creek, 49.2417°, −124.8000°: RBCM019936 (KF302820, NA, KF302970, KF3030480); (6) Canada, British Columbia: Vancouver Island, Lower Lost Shoe Creek, 49.0117°, −125.5894°: RBCM019803 (AY954931); RBCM019804 (KF302817, KF302889, KF302968, KF303046); (7) Canada, British Columbia: Vancouver Island, Veich Creek, 48.4333°, −123.6022°: RBCM019734 (KF302812, KF302886, KF302965, KF303043); Niagara Creek, 48.4939°, −123.5706°: RBCM019735 (KF302813, KF302887, KF302966, KF303044); Rithet Creek, 48.5839°, −123.7231°: RBCM019736 (KF302814); 48.5850°, −123.7231°: RBCM019737 (KF302815, KF302888, KF302967, KF303045).
Sorex navigator.—(8) Canada, British Columbia: Yoho River Valley, 51.4833°, −116.4667°: RBCM019354 (AY954932); Natural Bridge, 51.3833°, −116.5333°: RBCM019388 (AY954933); Leanchoil, 51.2167°, −116.5833°: RBCM019400 (AY954934); Sink Lake, 51.4500°, −116.3000°: RBCM019401 (AY954935); (9) Canada, British Columbia: Trojan Creek, 50.5522°, −120.9903°: RBCM020251 (KF302796, KF302872, KF302953, KF303030); Axe Creek, 50.5542°, −120.8675°: RBCM020252 (KF302797, KF302873, KF302954, KF303031); (10) Canada, British Columbia: Cayoosh Creek, 50.4986°, −122.2978°: RBCM019956 (AY954936, KF302864, KF302945, KF303022); RBCM019957 (AY954937, KF302865, KF302946, KF303023); (11) Canada, British Columbia: Callahan Creek, 50.1114°, −123.1114°: RBCM020315 (KF302798, KF302874, KF302955, KF303032); (12) Canada, British Columbia: Pinecone Burke Provincial Park, 49.4333°, −122.6667°: RBCM019823 (KF302785); RBCMMamm:019824 (KF302786); (13) United States, Washington: 48.9167°, −119.2833°: UWBM76620 (GU190749); (14) United States, Washington: 48.8333°, −117.1333°: UWBM77661 (EU856402); UWBM77663 (EU856404); (15) United States, Washington: 48.7333°, −120.6833°: UWBM76771 (EU856453); UWBM80443 (EU856421); (16) United States, Washington: 48.7000°, −121.1500°: UWBM74834 (EU856448); UWBM74895 (EU856449); (17) United States, Washington: 48.6833°, −120.6333°: UWBM80465 (GU190753); (18) United States, Washington: 48.6667°, −118.4333°: UWBM77659 (EU856400); (19) United States, Washington: 47.9667°, −123.2667°: UWBM75177 (EU856426); UWBM79799 (EU856454); 48.0000°, −123.6167°: UWBM79800 (EU856455); UWBM79807 (EU856456); 47.9833°, −123.6500°: UWBM79809 (EU856457); (20) United States, Washington: 47.8833°, −123.1500°: UWBM78827 (EU856410); (21) United States, Washington: 47.8333°, −120.6500°: UWBM75138 (EU856425); (22) United States, Washington: 47.7667°, −121.0500°: UWBM74821 (EU856447); (23) United States, Washington: 47.7000°, −123.9833°: UWBM80504 (EU856432); (24) United States, Washington: 47.5833°, −123.2500°: UWBM75686 (EU856443); 47.6000°, −123.0667°: UWBM80510 (GU190754); (25) United States, Washington: 46.3833°, −121.5833°: UWBM75841 (EU856444); 46.3000°, −121.7333°: UWBM80487 (EU856429); 46.2833°, −121.6500°: UWBM80489 (EU856431); (26) United States, Washington: 46.0833°, −117.8500°: UWBM77668 (EU856436); 46.1000°, −117.8500°: UWBM77669 (EU856437); (27) United States, Oregon: 45.2833°, −121.7500°: UWBM80401 (EU856415); 45.2000°, −121.6000°: UWBM80675 (EU856424); (28) United States, Oregon: 44.9500°, −118.2667°: UWBM77664 (GU190750); (29) United States, Oregon: Jefferson Co., Canyon Creek, 44.5104°, −121.7202°: MSB261725 (KF302780, KF302859, KF302940, NA); (30) United States, Oregon: Linn Co., Hackelman Creek, 44.3995°, −122.0892°: MSB261736 (KF302781, KF302860, KF302941, KF303018); (31) United States, Oregon: 44.0333°, −121.7500°: UWBM78799 (EU856407); (32) United States, Oregon: 42.7000°, −122.3333°: UWBM80595 (GU190755); 42.6667°, −122.3167°: UWBM80597 (GU190756); (33) United States, Alaska: Copper River Basin: MSB195275 (KF302774, KF302851, KF302932, KF303010); (34) United States, Alaska: Yakutat Quad: UAM76526 (KF302806, KF302881, KF302960, KF303038); (35) United States, Alaska: Circle, 65.2275°, −144.5003°: UAM34596 (KF302799); UAM34597 (KF302800); UAM34599 (KF302801, KF302876, NA, KF303034); UAM34600 (KF302802); UAM34601 (KF302803, KF302877, NA, NA); (36) United States, Alaska: Healy Quad, Moose Creek, 63.5000°, −150.5833°: UAM76302 (AF238034, KF302880, KF302959, KF303037); (37) United States, Alaska: Farm Island, Stikine River, 56.6333°, −132.4167°: UAM20787 (AF238033); Petersburg Quad, South of Thomas Bay, 57.0083°, −132.9833°: UAM24577 (NA, KF302956, KF302956, KF303033); (38) United States, Alaska: Wrangell Island, Pat's Lake, 56.3491°, −132.3383°: MSB195374 (KF302775, KF302852, KF302933, KF303011); Fools Creek, 56.2697°, −132.0706°: UAM49564 (KF302805, KF302879, KF302958, KF303036); (39) United States, Alaska: Bradfield Canal, 56.0269°, −130.0706°: UAM35350 (KF302804, KF302878, KF302957, KF303035); (40) Canada, Yukon: North Fork Klondike River, 64.0259°, −138.5789°: MSB145600 (KF302771, KF302848, KF302929, KF303008); (41) Canada, Yukon: Fox Creek, 61.1000°, −135.2934°: MSB144181 (KF302770, KF302847, KF302928, KF303007); (42) Canada, Yukon: Snafu Creek, 60.1446°, −133.8253°: MSB221077 (KF302776, KF302853, KF302934, KF303012); (43) Canada, Yukon: East of Watson Lake, 60.0651°, −128.6183°: MSB234567 (NA, NA, KF302939, KF303017); (44) Canada, British Columbia: Bulkley River, 55.2750°, −126.4617°: RBCM020007 (KF302787, KF302866, KF302947, KF303024); Twain Creek, 54.5433°, −125.9219°: RBCM020008 (KF302788, KF302867, KF302948, KF303025); 54.5781°, −125.9542°: RBCM020010 (KF302789, KF302868, KF302949, KF303026); RBCM020011 (KF302790, KF302869, KF302950, KF303027); (45) Canada, British Columbia: Morice River, 54.2072°, −127.1653°: RBCM020012 (KF302791); RBCM020017 (KF302795, KF302871, KF302952, KF303029); (46) Canada, British Columbia: Cluculz Creek, 53.7850°, −123.6931°: RBCM020013 (KF302792); 53.7756°, −123.6939°: RBCM020015 (KF302793); 53.7850°, −123.6931°: RBCM020016 (KF302794, KF302870, KF302951, KF303028); (47) United States, Idaho: Lemhi Co., Agency Creek, 44.9508°, −113.6134°: MSB227591 (NA, KF302856, KF302937, KF303015); 44.9759°, −113.5165°: MSB227766 (KF302779, NA, KF302938, KF303016); (48) United States, Idaho: Rock Creek, 42.2000°, −114.2823°: MSB225608 (KF302778, KF302855, KF302936, KF303014); (49) United States, Montana: Price Creek, 44.5659°, −112.1417°: MSB156230 (KF302772, KF302849, KF302930, NA); (50) United States, Montana: Hyalite Creek, 45.5541°, −111.0396°: MSB56695 (KF302782, KF302861, KF302942, KF303019); MSB56696 (KF302783, KF302862, KF302943, KF303020); (51) United States, Wyoming: Beartooth Hwy.: MSB123338 (KF302767, KF302843, KF302924, KF303003); Sunlight Creek: MSB123590 (KF302764, KF302844, KF302925, KF303004); (52) United States, Wyoming: Bighorn Mtns., Fool Creek, 44.8001°, −107.5727°: DNMS11044 (NA, JN889024, JN889244, JN889354); (53) United States, Wyoming: Squaw Creek, 43.1210°, −110.9255°: MSB155137 (JN889612, JN889025, JN889245, JN889355); (54) United States, Nevada: Bull Run Mtns., Rocky Gultch, 41.7773°, −115.9872°: MSB225575 (KF302777, KF302854, KF302935, KF303013); (55) United States, Utah: Zion National Park: MSB122504 (KF302766, KF302842, KF302923, KF303002); (56) United States, New Mexico: Jemez Mtns., Upper Puerco Crossing, 36.0870°, −106.7480°: MSB157027 (KF302773, KF302850, KF302931, KF303009); Jemez Mtns., 35.8148°, −106.5242°: MSB87740 (KF302784, KF302863, KF302944, KF303021); (57) United States, Arizona: White Mtns., 33.9571°, −109.5183°: MSB124084 (KF302768, KF302845, KF302926, KF303005); 33.9603°, −109.5061°: MSB124093 (KF302769, KF302846, KF302927, KF303006); (additional samples excluded from genealogies) Canada, British Columbia: 49.3667°, −122.9333°: UWBM73732 (EU856442); United States, Oregon: 42.7167°, −118.6167°: UWBM80543 (EU856397); UWBM80574 (EU856398); 44.0333°, −121.6667°: UWBM78773 (EU856406); 44.0667°, −121.6333°: UWBM78757 (EU856405); UWBM76723 (EU856452); 45.2333°, −121.7500°: UWBM80397 (EU856414); 45.3167°, −121.7167°: UWBM80368 (EU856411); UWBM80390 (EU856412); UWBM80394 (EU856413); United States, Washington: 46.0833°, −117.8500°: UWBM76688 (EU856451); UWBM77665 (EU856435); UWBM77671 (EU856439); 46.1000°, −117.8667°: UWBM77667 (GU190751); UWBM77670 (EU856438); UWBM77672 (EU856440); 46.1000°, −117.8333°: UWBM75236 (GU190748); UWBM74781 (GU190746); UWBM74782 (EU856445); UWBM74783 (EU856446); 46.1167°, −117.8333°: UWBM76682 (EU856450); 46.2833°, −121.6500°: UWBM80488 (EU856430); 46.3000°, −121.7333°: UWBM80486 (EU856428); 47.6000°, −123.0667°: UWBM80511 (EU856433); 47.8333°, −120.6500°: UWBM75141 (GU190747); 47.9667°, −123.1167°: UWBM78820 (EU856408); UWBM78823 (EU856409); 48.0500°, −121.7167°: UWBM80622 (EU856423); 48.6833°, −120.6333°: UWBM80448 (GU190752); 48.7000°, −120.6667°: UWBM80406 (EU856416); UWBM80410 (EU856417); UWBM80433 (EU856419); UWBM80464 (EU856422); 48.7333°, −120.6833°: UWBM80420 (EU856418); UWBM80442 (EU856420); 48.7333°, −118.1333°: UWBM77673 (EU856441); 48.7833°, −117.0833°: UWBM77662 (EU856403); 48.8000°, −117.0500°: UWBM77658 (EU856399); 48.8333°, −117.1833°: UWBM77660 (EU856401).
Sorex palustris.—(58) Canada, Alberta: Calling Lake, 55.1958°, −113.6583°: RBCM019829 (AY954939, KF302892, KF302971, KF303049); RBCM019830 (AY954938, KF302893, KF302972, KF303050); RBCM019831 (AY954940, KF302894, KF302973, KF303051); RBCM019832 (AY954941, KF302895, KF302974, KF303052); RBCM019833 (AY954942, KF302896, KF302975, KF303053); RBCM019837 (KF302827); (59) Canada, Alberta: Lac la Biche, 55.1333°, −111.6611°: RBCM019834 (AY954943, KF302897, KF302976, KF303054); RBCM019835 (KF302825, KF302898, KF302977, NA); RBCM019836 (KF302826, KF302899, KF302978, KF303055); RBCM019838 (KF302828, KF302900, KF302979, KF303056); (60) Canada, Alberta: Ghostpine Creek, 52.0319°, −113.3889°: RBCM019839 (KF302829, KF302901, KF302980, NA); RBCM019840 (KF302830, KF302902, KF302981, NA); Red Deer River, 52.0625°, −114.0525°: RBCM019841 (KF302831, KF302903, KF302982, NA); (61) Canada, Alberta: Okotoks: ROM109832 (DQ788800); (62) Canada, Manitoba: Caddy Lake (DQ995274); (63) Canada, Manitoba: Piney (DQ995275).
Sorex albibarbis.—(64) Canada, Nova Scotia: Cape Breton (DQ995278); Cape Breton Islands (DQ788801); (65) Canada, Nova Scotia: Mount Martock (DQ995276); Waverly (DQ995277); Cole Harbor: ROM110278 (DQ788802); (66) Canada, New Brunswick: Edmundston (DQ995279; DQ995280); (67) Canada, Quebec: 50.7353°, −74.9707°: MSB229151 (KF302765, KF302904, KF302983, KF303057); (68) Canada, Newfoundland and Labrador: North West River, 53.5683°, −60.0800°: MSB229139 (KF302832, KF302905, KF302984, KF303058).
Sorex bendirii.—(69) Canada, British Columbia: Seymour River, 49.3500°, −123.0000°: RBCM019527 (AY954946); Delta, 49.1350°, −123.0083°: RBCM020945 (KF302835, KF302912, KF302991, NA); Squamish, 49.7369°, −123.1103°: RBCM020946 (KF302836, KF302913, NA, NA); Surrey, Hjorth Creek, 49.1911°, −122.7742°: RBCM021087 (NA, KF302914, KF302993, NA); (70) Canada, British Columbia: Sumas Mountain, Clayburn Creek, 49.0833°, −122.2167°: RBCM019434 (AY954944); Sumas Mountain, McKee Creek, 49.0667°, −122.2833°: RBCM019435 (AY954945); RBCM019436 (AY954947, KF302908, KF302987, NA); Aldergrove, 49.0839°, −122.5014°: RBCM019822 (NA, KF302909, KF302988, KF303061); Hemlock Valley, 49.3258°, −121.8625°: RBCM020943 (KF302833, KF302910, KF302989, NA); 49.3261°, −121.8611°: RBCM020944 (KF302834, KF302911, KF302990, NA); OPSEE Military Reserve, 49.0672°, −121.9308°: RBCM020947 (NA, NA, KF302992, NA); (71) United States, California: Southeast of Trinidad, 41.0040°, −124.0916°: MSB43606 (KF302762, KF302906, KF302985, KF303059); (72) United States, Oregon, Tillamook Co., 45.2669°, −123.8836°: UAM52163 (AF238031); (73) United States, Oregon: Lane Co., 44.2417°, −123.8250°: UAM52161 (AF238032).
Sorex cinereus.—Canada, British Columbia: Nass River Valley, 55.3319°, −128.9653°: UAM52330 (JN889490, JN888994, JN889214, JN889324).
Sorex fumeus.—Canada, Ontario: Northwest of Gilmour, 40.1568°, −79.2578°: MSB53189 (KF302841, KF302921, KF303000, KF303068); United States, Pennsylvania: Powdermill Nature Reserve, 47.7191°, −122.7185°: MSB53147 (JN889645, KF302920, KF302999, KF303067); United States, Massachusetts: Northwest of Whately Center: MSB47899 (KF302840, KF302919, KF302998, KF303066).
Sorex haydeni.—United States, South Dakota: Bluestem Game Production Area, 43.6700°, −98.1400°: UAM50345 (AY014939, JN889004, JN889224, JN889334).
Sorex monticolus.—United States, Colorado: Apishapa Creek, 37.3411°, −105.0006°: DNMS11006 (JN889534, JN889019, JN889239, JN889349); United States, Colorado: Commisary Creek, 38.5608°, −107.3247°: MSB89327 (KF302763, KF302922, KF303001, NA); United States, Montana: Price Creek, 44.5659°, −112.1417°: MSB156241 (JN889535, JN889021, JN889241, JN889351); United States, New Mexico: Jemez Mtns., 44.8008°, −77.6720°: MSB140765 (KF302837, KF302915, KF302994, KF303062); United States, Alaska: Palmer Creek Road, 60.8312°, −149.5359°: MSB145441 (JN889608, JN889020, JN889240, JN889350).
Sorex ornatus.—United States, California: Heart Bar Campground, 34.1586°, −116.7860°: MSB40419 (KF302839, KF302917, KF302996, KF303064).
Sorex pacificus.—United States, Oregon: Benton Co., Alsea Falls, 44.3225°, −123.4752°: MSB43569 (AF238020, KF302918, KF302997, KF303065).
Sorex rohweri.—United States, Washington: Jefferson Co., Hoh Valley, 47.8300°, −123.9800°: UWBM74992 (JN889615, JN889056, JN889276, JN889386).
Sorex sonomae.—United States, California: North of Trinidad, 41.1362°, −124.1548°: MSB43641 (AF238026, KF302907, KF302986, KF303060).
Sorex trowbridgii.—United States, Washington: Naval Sub Base, Bangor: MSB83485 (JN889622, JN889058, JN889278, JN889388).
Sorex vagrans.—Canada, British Columbia: Blueberry Creek, 49.2580°, −117.9395°: MSB156389 (KF302838, KF302916, KF302995, KF303063).
GenBank sequences originally cited in Demboski and Cook (2001): AF238020, AF238026, AF238031–AF238034, AY014939; O'Neill et al. (2005): AY954926–AY954947; Shafer and Stewart (2007): DQ788800–DQ788802; Mycroft et al. (2011): DQ995274–DQ995280; Himes and Kenagy (2010): EU856397–EU856426, EU856428–EU856433, EU856435–EU856457, GU190746–GU190756; Hope et al. (2012): JN889490, JN889534, JN889535, JN889608, JN889612, JN889615, JN889622, JN889645, JN888994, JN889004, JN889019–JN889021, JN889024, JN889025, JN889056, JN889058, JN889214, JN889224, JN889239–JN889241, JN889244, JN889245, JN889276, JN889278, JN889324, JN889334, JN889349–JN889351, JN889354, JN889355, JN889386, JN889388.
