OUP user menu

Phylogeography of the western jumping mouse (Zapus princeps) detects deep and persistent allopatry with expansion

Jason L. Malaney, Chris J. Conroy, Lena A. Moffitt, Harmony D. Spoonhunter, James L. Patton, Joseph A. Cook
DOI: http://dx.doi.org/10.1644/12-MAMM-A-006.1 1016-1029 First published online: 15 October 2013

Abstract

Understanding how diversity is partitioned across the landscape provides perspectives on the environmental processes that have influenced the evolutionary history of organisms. We analyzed spatial demography, historical biogeography, and niche divergence of the western jumping mouse (Zapus princeps) using molecular sequence data from mitochondrial and nuclear DNA recovered from 7 of the 11 subspecies in western North America. Phylogeographic structure within Z. princeps was partitioned across 5 clades (Boreal, Northern Sierra, Southern Rockies, Southern Sierra, and Uinta). Two lineages detected in the Sierra Nevada of California (Northern Sierra and Southern Sierra) were more closely allied to Z. trinotatus than to other lineages of Z. princeps and species distribution models mirror these phylogenetic signatures by detecting wide overlap in niches for Sierran jumping mice and Z. trinotatus as compared to other Z. princeps. Four southern lineages are deeply divergent and limited to disjunct mesic and montane habitats within the xeric southwestern United States, whereas the 5th lineage is widespread, extending from Wyoming to Alaska and reflecting expansion northward following deglaciation, a common pattern in boreal mammals.

Key words
  • genetic diversity
  • jumping mice
  • Rocky Mountains
  • Sierra Nevada
  • species distribution model
  • Zapus princeps

Western North America has a diverse biota that is the product of complex evolutionary and environmental processes (Lomolino et al. 2006). Significant intraspecific genetic variation in mammals in the region is hypothesized to have been shaped by extreme topographic heterogeneity and repeated glaciations during the Quaternary (Riddle and Hafner 2006). Molecular investigations of western mammals have provided new views of diversification, occasionally revealing unexpected genetic architectures (Riddle et al. 2000; Matocq 2002b; Alvarez-Castañeda and Patton 2004; Galbreath et al. 2010). Documentation of geographic molecular variation in organisms, when combined with assessments of demography, historical biogeography, and niche variation, can provide insight into key questions related to climate change, post-Pleistocene colonization, habitat fragmentation, and possible future response to changing environments (Avise 2000; Riddle and Hafner 2007). We assess phylogeographic structure in the western jumping mouse (Zapus princeps) to further explore the biogeographic history of montane regions in western North America.

The western jumping mouse is an inhabitant of mesic and montane habitats ranging from New Mexico and central California northward through most of western North America (Fig. 1) to southeastern Alaska and southern Yukon, Canada (Krutzsch 1954; Hafner et al. 1981; Hall 1981; Jones 1981). Because much of this widespread range was glaciated during the late Quaternary, the patterns and levels of connectivity of the mesic environments these mice inhabit have changed, potentially leaving an imprint on molecular variation in the species (Hewitt 1996, 2000; Waltari and Guralnick 2009). The dynamic geologic history, variable topography, and patchy mesic environments of western North America provide a series of evolutionary experiments because multiple mountain ranges may represent replicated isolation events. We explore historical-biogeographic questions related to the effects of Pleistocene fragmentation on mesic-associated biota by testing 2 primary ideas using molecular variation and niche modeling in jumping mice. First, we explore lineage divergence related to glacial cycling by asking whether jumping mice lineages were in wide contact during glaciations or if they remained geographically isolated, similar to their current distribution. Second, given that their contemporary northern range was blanketed by ice until the end of the Pleistocene, we test the scenario that southern lineages of Z. princeps reflect signatures of genetic structure that are deeper (due to persistence) than northern lineage(s), which presumably expanded (ephemeral) from the south following the latest Pleistocene deglaciation.

Fig. 1

Distribution and range limits of Zapus princeps (stipple) and Z. trinotatus (hash) modified from Hall (1981). Currently recognized subspecies: 1) Z. p. chrysogenys, 2) Z. p. cinereus, 3) Z. p, curtatus, 4) Z. p. idahoensis, 5) Z. p. kootenayensis, 6) Z. p. minor, 7) Z. p. oregonus, 8) Z. p. pacificus, 9) Z. p. princeps, 10) Z. psaltator, 11)Z.p. utahensis, a) Z. t. trinotatus, b) Z. t. montanus, c) Z. t. eureka, and d) Z. orarius.

Two mutually exclusive and competing hypotheses, related to lineage divergence and historical biogeography, are plausible. The admixture scenario, with populations of jumping mice isolated during interglacials on montane islands (e.g., the contemporary condition) but in wide contact during glacial advances, predicts that genetic signatures should reflect high levels of exchange. Previous studies in western montane environments have identified different forms and levels of mixing (admixture, introgression, and hybridization) in grasshoppers (Knowles 2001), birds (Spellman and Klicka 2007; Mettler and Spellman 2009), pikas (Galbreath et al. 2009, 2010), and rodents (Good and Sullivan 2001; Spaeth et al. 2009). Alternatively, jumping mice lineages may have remained independent and evolved in situ due to sustained geographic isolation or persistent allopatry. Evidence for the persistent allopatry scenario would include geographically structured lineages. Signatures of genetic divergence and geographically structured variation in montane organisms have been recorded for alpine stonecrops (DeChaine and Martin 2005b), foxtail pine (Eckert et al. 2008), kittentails (Marlowe and Hufford 2008), alpine butterflies (DeChaine and Martin 2005a), birds (Spellman et al. 2007), and various mammals (Sullivan et al. 2000; Demboski and Cook 2001; Hornsby and Matocq 2011).

North American environments have not been static and organisms have presumably tracked climate change (Hewitt 2000), with many species dispersing northward following glacial retreat (Anderson and Boms 1994; Lessa et al. 2003). In Canada and Alaska, northern populations of Z. princeps are likely the result of post-Pleistocene colonization of deglaciated terrains, as observed in other mammals (Arbogast 1999; Conroy and Cook 2000b; Runck and Cook 2005). Southwestern peripheral populations of jumping mice, in contrast, are comparatively more fragmented and largely restricted to dense alder (Alnus spp.), willow (Salix spp.), and aspen (Populus spp.) stands typically associated with riparian systems and high-elevation mesic habitats (Quimby 1951; Krutzsch 1954; Hart et al. 2004; Frey and Malaney 2009). Because these habitats are limited in extent and isolated within predominantly xeric environments, populations of montane mesic-associated mammals are hypothesized to exhibit a deeper divergence than their northern counterparts. For example, genetic breaks across arid barriers such as the Columbia and Wyoming basins (Nielson et al. 2001; Carstens et al. 2005; DeChaine and Martin 2005b; Carstens and Richards 2007) suggest that isolation during warm interglacial periods may contribute to allopatric divergence in mammals (Demboski and Cook 2001; Demboski and Sullivan 2003; Good et al. 2008; Galbreath et al. 2010; Homsby and Matocq 2011). Persistent allopatric divergence is expected to result in deeper genetic distance and lineage cohesion, suggesting limited or no mixing through multiple glacial cycles (Nielson et al. 2001, 2006; DeChaine and Martin 2006).

We begin by assessing how past events may have influenced the phylogenetic signature, demography, and historical biogeography of the western jumping mouse. We test if the phylogenetic signal is the result of admixture or persistent allopatry using coalescent simulations. Next, we examine the existing taxonomic framework for Z. princeps, originally based on morphological features (Krutzsch 1954), and then assess environmental variation with species distribution models and test for niche and range overlap among paraphyletic clades. Finally, we propose a set of alternative biogeographic hypotheses based on allopatric modes of speciation that may account for the observed phylogenetic signal.

Materials and Methods

Extraction and sequencing of DNA.— We included specimens from the Museum of Southwestern Biology (MSB—University of New Mexico, Albuquerque), the Museum of Vertebrate Zoology (MVZ—University of California, Berkeley), and the University of Alaska Museum of the North (UAM—University of Alaska, Fairbanks; Supporting Information S1, DOI: 10.1664/12-MAMM-A-006.S1). Specimens represent 7 of the 11 subspecies (Z. p. idahoensis, Z. p. minor, Z. p. oregonus, Z. p. pacificus, Z. p. princeps, Z. p. saltator, and Z. p. utahensis) recognized by Krutzsch (1954) and Hall (1981) and range from southeastern Alaska and Yukon southward throughout most western states to the species' southem limits in Califomia and New Mexico (Fig. 1). Detailed information related to voucher specimens is available through the Arctos database (Arctos Collection Management Information System 2010). The mitochondrial DNA (mtDNA) cytochrome-b (Cytb, 1,140 base pairs [bp]) gene was obtained for 91 specimens (see Supporting Information S1) representing 46 locations across the range of Z. princeps. In addition, samples of Z. trinotatus (4), Z. hudsonius (6), and Napaeozapus insignis (1) were included. Independent evolutionary perspectives were gained by sequencing 2 nuclear DNA markers to test major lineage breaks identified by mtDNA; glucocerebrosidase (GBA, 347 bp) and myosin heavy chain 2 (MYH2, 267 bp).

Total genomic DNA was extracted from frozen or ethanolpreserved tissues (heart or liver). Amplification of the Cytbgene was conducted with primers L14724 and H15915 (Irwin et al. 1991) or with a combination of MVZ05–MVZ16 and MVZ127–MVZ108 (Smith and Patton 1993; Leite and Patton 2002) using protocols previously established (Lessa and Cook 1998; Halanych et al. 1999; Patton et al. 2008). Polymerase chain reaction products were sequenced using BigDye Terminator Cycle Sequencing Ready Reaction mix (version 3.1; Applied Biosystems, Foster City, California) with combinations of amplification primers. Heavy and light strands were sequenced in both directions using an Applied Biosystems 3100 automated DNA sequencer (Applied Biosystems) in the Molecular Biology Facility, Biology Department, University of New Mexico, Albuquerque, or at the Museum of Vertebrate Zoology, University of California, Berkeley. Sequences were analyzed using Sequencher (version 4.9; Gene Codes Corporation, Inc., Ann Arbor, Michigan) with a reference sequence from GenBank (Z. trinotatus, AF1 19262).

We sequenced 2 fragments of nuclear introns (GBA and MYH2) for several specimens using published primers and protocols (Lyons et al. 1997). We sequenced a minimum of 2 randomly selected individuals from mitochondrial clades of Z. princeps (30 GBA and 14 MYH2), Z. trinotatus (2), and Z. hudsonius (2) with 1 N. insignis used to root phylogenies. All individuals were sequenced in both directions. Heterozygous positions were identified and polymorphic alleles were assessed using Phase (version 2.1—Stephens et al. 2001; Stephens and Scheet 2005) in DnaSP (version 5.0—Rozas et al. 2003; Librado and Rozas 2009) with haplotypes inferred from multiallelic loci using a Bayesian framework with 0.90 cutoff and 10,000 iterations. Unresolved haplotypes were coded as missing data.

Alignments were completed using default parameters and algorithms of Clustal X (Larkin et al. 2007) in the program MEGA (version 5.05—Tamura et al. 2011). Contigs for all genes are available in GenBank (Supporting Information S1).

Diversity measures.—Molecular diversity and demographic estimates from each marker were determined for putative Z. princeps (Table 1). Neutrality and population equilibrium were assessed via Tajima's D and Fu's FS tests and 10,000 coalescent simulations to assess significance. The mtDNA data set was partitioned into clades to assess demographic change. Net sequence divergence (da) was calculated between the observed mtDNA clades (Nei 1987). We calculated segregating sites (S), haplotypes (h), haplotype diversity (ƞ), and nucleotide diversity (Л) for each marker with DnaSP (Rozas et al. 2003; Librado and Rozas 2009).

View this table:
Table 1

Molecular diversity indexes of mitochondrial DNA (cytochrome-b [Cytb]) and nuclear DNA (glucocerebrosidase [GBA] and myosin heavy chain 2 [MYH2]) for North American jumping mice calculated in DnaSP; n = number of individuals sampled, S = polymorphic sites, h = number of haplotypes, ƞ = haplotype diversity and standard deviation (SD), л = nucleotide diversity (per site), and neutrality estimates. Tajima's D and Fu's FS were calculated with coalescent simulations (10,000 replicates) with values that correspond to calculations based on the Waterson estimator theta (ϴw). Asterisks represent significance: * < 0.05; ** < 0.01. Napaeozapus was excluded from Cytb analyses.

nShƞ (SD)лDFs
Cytb gene101157410.933 (0.018)0.0931.3036.981
Boreal3275140.986 (0.013)0.010-1.710*-11.732**
Northern Sierra51151.000 (0.126)0.005
Southem Rockies143070.846 (0.061)0.0110.9924.720
Southem Sierra361550.546 (0.094)0.005-1.626*-5.789**
Uinta4330.833 (0.222)0.002
Z. trinotatus4203
Z. hudsonius6824
GBA gene32 (64)22180.8950.010-1.030-5.893*
MYH2 gene25 (50)16170.9380.014-0.305-9.214*

Phylogenetic and phylogeographic analyses.—Our aligned mitochondrial and nuclear data were processed via MrModeltest (Nylander 2004) performing hierarchical likelihood-ratio tests and calculating Akaike information criterion (AIC); both measures agreed for all genes. General time reversal (Tavare 1986) plus gamma (1.6506) plus proportion (0.4875) of invariant sites model (GTR+Γ+I, log-likelihood = -5,760.9790, K = 10, AIC = 11,541.9580) was selected as the most appropriate evolutionary model for the mitochondrial marker and subsequently used in Bayesian inference and maximum-likelihood analyses. The Kimura (1980) model (K80, log-likelihood = -661.0489, K = 1, AIC = 1,324.0778) was selected for the GBA gene. The Hasegawa, Kishino, and Yano (Hasegawa et al. 1985) plus proportion (0.9085) invariant sites model (HKY+I, log-likelihood = 436.1803, K = 5, AIC = 882.3607) was identified for the MYH2 gene.

Bayesian reconstruction was performed using MRBAYES (version 3.1.2—Huelsenbeck and Ronquist 2001; Lakner et al. 2008) beginning with random trees and Markov chain sampled every 1,000th tree for 2 million generations and 4 chains run simultaneously with temperature set to 0.20 for 3 chains and 1 cold chain. Three replicate runs were completed to confirm consistency and each marker was run with distinct priors set from MrModeltest output. Chain stationarity was assessed by inspecting the standard deviation of split frequencies consistently below 0.05 and confirmed complete via the graphical output from the initial 50,000 generations with 0.20 of each replicate discarded as bum-in (Huelsenbeck and Imennov 2002). Nodal strength (posterior probability) was identified in the consensus of the residual trees and the midpoint rooted majority-rule consensus tree was visualized in the program FigTree version 1.3.1 (FigTree 2010; Fig. 2).

Fig. 2

Phylogram of Bayesian inference majority rules consensus tree produced using GTR+Γ+I model and samples of Zapus princeps and other taxa for the mitochondrial cytochrome-b gene. Asterisks at nodes correspond to posterior probabilities (> 0.95, PP) from 50,000 post-burn-in trees and 1,000 nonparametric bootstraps (> 0.90, maximum likelihood [ML]).

Maximum-likelihood optimality criteria were used for phylogenetic reconstruction using garli (version 2.0 parallel—Zwickl 2006). We considered all characters as unordered with 4 possible states (A, C, G, T) with heuristic searches. Distinct models of evolution were applied to each marker with discrete base frequencies and rate categories for each from MrModeltest. Tree-bisection-reconnection branch swapping was employed with 100 random stepwise additions. Three runs were conducted to ensure consistency and nonparametric bootstrap support (Felsenstein 1985) was evaluated with 1,000 pseudoreplicates (Fig. 2).

The median-joining statistical parsimony algorithm (Bandelt et al. 1999) in the program network (version 4.2; Fluxus Engineering, Suffolk, United Kingdom) was employed to produce a haplotype network (Fig. 3) for each marker given that intraspecific phylogenetic methods may fail (Posada and Crandall 2002). This algorithm calculates the similarity between haplotypes into a network where the combined probability is > 95% (Templeton et al. 1992).

Fig. 3

Median-joining statistical parsimony networks for jumping mice mitochondrial DNA (cytochrome-b [Cytb]) with respect to geography and nuclear DNA (glucocerebrosidase [GBA] and myosin heavy chain 2 [MYH2]—inset). Individual tick marks represent 1 polymorphic site or mutation (step) and squares represent missing or ancestral haplotypes. Note multiple haplotypes per locality in southern lineages reflecting deep phylogenetic history but haplotypes shared among geographic locations in the most northern populations of the Boreal clade reflecting recent demographic history.

Taxonomic evaluation.—We compared the morphological taxonomic classification of Z. princeps (Krutzsch 1954) against the mitochondrial framework using the Shimodaira—Hasegawa test (Shimodaira and Hasegawa 1999). A maximum-likelihood tree constrained to reflect monophyly of Z. princeps, Z. trinotatus, and Z. hudsonius was compared to the unconstrained best maximum-likelihood tree using PAUP* (version 4.0b 10—Swofford 2003), GTR+Γ+I model of nucleotide substitution, and 10,000 resampling of estimated log-likelihoods bootstrap replicates (Hasegawa and Kishino 1994).

Coalescent simulations.—We conducted coalescent simulations using the parametric bootstrap method (Goldman et al. 2000). Our aim was to test altemative hypotheses of admixture versus persistent allopatry using a likelihood-based or frequentist approach (Knowles and Maddison 2002; Hickerson et al. 2010) with coalescent simulations in the program MESQUITE (Maddison and Maddison 2009). Alternative hypothesized phylogenies (admixture versus persistent allopatry) were simulated for 1,000 replicates to produce gene matrices using ancestral NE(f) (183,779) from our estimate of the parameter ϴW (calculated in DnaSP), branches scaled to branchhig pattern, 1-year generation time (Brown 1970; Cranford 1983), and ancestral divergence estimates from zapodid fossil records (Kurtén and Anderson 1980; Hafner 1993; Ruez and Bell 2004). To assess the validity of each model and determine overlap between models, we used upper and lower confidence intervals (90%) of ϴW (Knowles and Carstens 2007) on independent runs.

From the 1,000 gene matrices, we constructed genealogies in PAUP* using heuristic parsimony searches with 10 random addition replicates and tree-bisection-reconnection branch swapping with max-trees set to 100, and produced a majority-rule consensus tree. Next we calculated the discord between the reconstructed gene tree and the assignment of individuals into separate lineages using s of Slatkin and Maddison (1989) as implemented in MESQUITE. Finally, assessment of our 2-tailed test was considered significant if the empirical data occur outside of the 90% confidence interval of the simulation data.

Our ϴ estimates (effective population size scaled to the neutral mutation rate) were calibrated to recent fossil dates using the equation ϴ = 4NE(f)µ, assuming µ = 3.14 substitutions per million years as calculated by Malaney et al. (2012) and scaled branch widths of ϴ were used for the competitive hypotheses.

Species distribution modeling.—Climate variables such as temperature and precipitation are known to affect the metabolic rates of jumping mice (Cranford 1975) and can provide insight into the spatial distribution of environmental characteristics for monophyletic lineages. We obtained bioclimatic variables from 2.5-min (4-km) resolution coverages from the Worldclim database (Humans et al. 2005; Worldclim 2005). Torpor in jumping mice is impacted by elevation (Cranford 1978; Muchlinski and Rybak 1978), so we included this coverage.

We followed modeling procedures from previous studies (Waltari et al. 2007; Waltari and Guralnick 2009) by clipping the coverages to the study area (western North America) and reducing the data set (Rissler and Apodaca 2007) to the 12 most biologically meaningful and uncorrelated coverages (Bio1—Annual Mean Temperature, Bio2—Mean Diurnal Range, Bio3—Isothermality, Bio7—Temperature Annual Range, Bio8—Mean Temperature of Wettest Quarter, Bio9— Mean Temperature of Driest Quarter, Bio 15—Precipitation Seasonality, Bio 16—^Precipitation of Wettest Quarter, Bio 17— Precipitation of Driest Quarter, Bio 18—Precipitation of Warmest Quarter, Bio 19—Precipitation of Coldest Quarter, and elevation). Localities for Z. princeps and Z. trinotatus were downloaded from MaNIS (Manis 2010). Localities with > 0.5 km2 uncertainty were discarded and several records were georeferenced (JLM) using BioGeomancer (BioGeomancer 2006; Guralnick et al. 2006). To account for sampling biases (Reddy and Davalos 2003), which may result in model overfitting and subjective outcome, we removed redundant records and then down-sampled localities to 10-km intervals by removing intervening records. To test if the Sierra Nevada lineages reflect phylogenetic signal in niche occupancy, we also partitioned localities of Z. princeps (see below). We constructed species distribution models using the default settings in the program Maxent (version 3.3.3a—Elith et al. 2006; Phillips et al. 2006) and ran 20 replicates with randomized 20th percentile training presence and depicted results using the pointwise bootstrap mean.

To identify overlap between taxonomic divisions and establish whether Sierra Nevada clades occupy analogous environments to Z. trinotatus or other Z. princeps, we completed a series of niche tests. Localities were partitioned into 3 groups reflecting DNA phylogenetic signal: Z. trinotatus, Sierran (Southem Sierra and Northern Sierra) clades, and remaining Z. princeps (Boreal, Southem Rockies, and Uinta) clades. First, we calculated the proportion of pixels (km2) where overlap between suitable niches occur using ArcGIS (version 10.0— ESRI 2010). Threshold values were determined from the more conservative “last sample included” criterion and were 0.21 for Z. princeps, 0.24 for Z. trinotatus, and 0.28 for Sierran. Quantifying niche overlap was accomplished with 3 metrics using ENMTools (version 3.1—Warren et al. 2010); Schoener's D (Schoener 1968), Warren's I (Warren et al. 2008), and relative ranks (RRWarren and Seifert 2011). Each measure identifies pairwise niche overlap values between 0 (none) and 1 (full). We calculated 100 pseudoreplicates of niche models (Warren et al. 2008, 2010) and corresponding measures of niche overlap among lineages following methods of Pyron and Burbrink (2009). This conforms to a 1-tailed test to identify if niche models are significantly different from a null distribution by randomly assigning lineage membership to the occurrence points for any 2 lineages. All species distribution models have basic assumptions including niche conservatism (Wiens and Graham 2005), whether coverages (environmental data) are adequate to generate predictions of a species' distribution (Kozak et al. 2008; McCormack et al. 2010), and adequate occurrence points to encapsulate the range of environmental conditions in the species niche (Pearson et al. 2007).

Results

Molecular diversity.—The Cytb gene was obtained for 91 specimens of Z. princeps, 4 Z. trinotatus, 6 Z. hudsonius, and 1 N. insignis (Supporting Information SI). Excluding N. insignis, there were 157 segregating sites (Table 1). Nucleotide composition (28.4% adenine, 32.0% thymine, 26.9% cytosine, and 12.7% guanine), transition : trans version ratio (R = 4.73), and codon position changes (6—1st position, 0— 2nd position, and 39—3rd position) among Zapus lineages were consistent with other measures of genuine Cytb gene for mammals (Irwin et al. 1991).

Twenty-seven (GBA) and 20 (MYH2) randomly selected individuals of Z. princeps from the 5 mitochondrial clades were sequenced for the nuclear introns (Lyons et al. 1997) plus 2 Z. hudsonius, 2 Z. trinotatus, and 1 N. insignis. For the GBA marker nucleotide composition was 21.6% adenine, 24.1% thymine, 28.3% cytosine, and 26.0% guanine. There were 22 segregating sites, 1 site with more than 2 variants, transition : transversion ratio of 2.26, plus a 2-bp insertion—deletion between Zapus and Napaeozapus. The MYH2 marker had 20.3% adenine, 25.2% thymine, 28.1% cytosine, and 26.4% guanine, 16 segregating sites, and 2 multivariant sites. Transition : transversion ratio was 2.11 with 9 positions represented by insertion—deletions.

Phylogenetic and phylogeographic divergence.— Phylogenetic reconstructions using maximum-likelihood and Bayesian (posterior probability) techniques and the parsimony network indicated congruent topologies (Figs. 2 and 3) for all genes. Seven clades were recovered including discrete Z. trinotatus and Z. hudsonius, plus 5 putative Z. princeps clades (Boreal, Northern Sierra, Southern Rockies, Southem Sierra, and Uinta [Figs. 2 and 3]). Each of the 7 clades identified in the phylogram (Fig. 2) reflects high bootstrap support (maximum likelihood) and posterior probabilities for major nodes. Both of the nuclear perspectives reflect deep nodes but provided less resolution near the tips. Current taxonomy of Z. princeps (Holden and Musser 2005) showed a polyphyletic relationship with respect to Z. trinotatus based on both the mitochondrial and nuclear data, with Northern Sierra and Southem Sierra clades more closely related to Z. trinotatus than to other Z. princeps (Boreal, Southem Rockies, and Uinta [Figs. 2 and 3]).

There were 42 Cytb haplotypes in our data set (Hd = 0.993; Fig. 3) with 14 haplotypes representing the Boreal clade, 5 in the Northern Sierra, 7 in the Southem Rockies, 5 in the Southem Sierra, 3 in the Uinta, and 3 in Z. trinotatus. Five other haplotypes including 4 in Z. hudsonius and 1 in N. insignis were identified (haplotypes 38—42 not displayed in Fig. 3). No widespread haplotype sharing was documented (Figs. 2 and 3) but a few haplotypes were shared between adjacent sampling localities within clades. Further, no haplotypes were colocated at a single locality from distinct clades. Several polymorphic sites define each lineage of jumping mice, with highest geographic structure across the 4 southem lineages (Northern Sierra, Southem Rockies, Southern Sierra, and Uinta [Fig. 3]). For example, > 100 steps separated the Uinta from the Southem Rockies and Boreal clades, whereas at least 47 steps segregated haplotypes of the Boreal and Southem Sierra clades. The Northern Sierra clade and Z. trinotatus are segregated by < 60 mutational steps. Between Z. princeps (Uinta) and Z. trinotatus there are 136 steps, whereas there are at least 147 mutational steps between Z. trinotatus and Z. hudsonius. Almost 300 mutational steps separated Z. hudsonius (haplotype 2; Z. h. luteus) and the Uinta Z. princeps (haplotype 51). Haplotype structure for the independent nuclear genes was less pronounced (fewer segregating sites) but consistent with the mtDNA signature.

Neutrality tests for mitochondrial lineages (Tajima's D and Fu's FS, Table 1) for the Boreal and Southem Sierra clades were significantly negative, suggesting deviation from mutation-drift equilibrium and may suggest population expansion (Excoffier et al. 2009). Results of neutrality tests for the Southem Rockies clade were not significant and small sample sizes from the Northern Sierra and Uinta clades precluded population-level analyses.

Intraclade diversity.—Patterns of intraclade genetic diversity differed across clades (Table 1). Within the Boreal and Southem Sierra clades, haplotypes showed a “starlike” phylogeny (Figs. 2 and 3) with relatively few mutational steps between haplotypes. In contrast, the Southem Rockies clade showed higher structure and a greater number of mutational steps. For example, within the Southem Rockies there were 16 steps separating haplotype 16 (Jackson County, Colorado) and haplotype 21 (Santa Fe County, New Mexico) over a geographic distance of 500 km. In the Boreal clade there were only 8 steps between haplotypes 11 (Yellowstone National Park) and 4 (Unuk River, Alaska) over a geographic distance of > 2,000 km (Fig. 3). The Uinta clade was represented by 4 specimens from the same locality (Strawberry Reservoir) and 3 haplotypes with 5 mutations. The 5 haplotypes from the Northern Sierra clade were distributed among 3 localities segregated by 4 mutational steps over 250 km. The Southem Sierra clade was represented by 11 sampling localities and 5 haplotypes within 200 km.

Taxonomic evaluation and coalescent simulations.—A tree constrained for monophyly of Z. princeps was significantly worse (P < 0.001) than the unconstrained maximum likelihood topology that showed paraphyly of these lineages (maximum likelihoodbest = -4,635.7099; maximum likelihoodconstraints = 8,266.8500 [Shimodaira and Hasegawa 1999]). Results of the parametric bootstrap test of altemative hypotheses confirmed the persistent allopatry hypothesis as the best match to the empirical data (Fig. 4). The empirical data had an s-value of 13 that was significantly different (2-tailed) from the admixture hypothesis.

Fig. 4

Altemative demographic hypotheses for jumping mice clades using the parametric bootstrap test (coalescent simulations) of divergence patterns with 90% confidence intervals. Hypotheses of persistent allopatry (PA, gray) versus admixture (AM, white) with the arrow highlighting the empirical tree value for s of Slatkin and Maddison s (1989) (s = 13). Numbers correspond with geographic ranges and clades in PA (1 = Z. trinotatus, 2 = Northern Sierra, 3 = Southem Sierra, 4 = Boreal, 5 = Uinta, 6 = Southern Rockies) but clades are represented geographically (mixed) in AM.

Species distribution models.—Modeling procedures had high area under the curve scores (> 0.95 for each model). Models were based on 170 localities for Z. trinotatus and 499 putative localities for Z. princeps; of these, 66 Sierra Nevada localities (Z. p. pacificus) were partitioned from other lineages of Z. princeps. Niche overlap and range overlap (threshold value > 0.20) was higher for the Sierran jumping mice and Z. trinotatus than the Sierran jumping mice and other Z. princeps (Fig. 5). All pairwise comparisons using pseudoreplicates reflected values significantly different than predicted from random among lineages.

Fig. 5

Phylogenetically informed species distributions models using Maxent based on pointwise mean logistic bootstrap prediction from 20 replicates and minimum training presence threshold rule (> 0.20) from 12 environmental variables for Zapus princeps, Z. trinotatus, and Sierran jumping mice. Geographic overlap was calculated in square kilometers for each pair and indexes (D, I, and RR—Warren et al. 2008, 2010) with significance (* < 0.05; ** < 0.001) via 100 pseudoreplicates.

Discussion

We documented a history of demographic and range expansion for the Boreal clade that contrasts with long-term persistent southem clades (Northern Sierra, Southem Rockies, Southem Sierra, and Uinta). This spatial and temporal contrast reflects distinctive demographic processes across the distribution of jumping mice of the west. The phylogeny and wide niche overlap between Sierran jumping mice and Z. trinotatus are inconsistent with current taxonomy that identifies the former as Z. princeps (Hall 1981; Krutzsch 1954; Holden and Musser 2005). Finally, the coalescent-based tests provide clarity on the historical biogeography of jumping mice over multiple glacial cycles.

Spatial demography.—For species with broad latitudinal distributions in North America, a common partem of deep southern (persistent) and shallow northern (ephemeral) stmcture has been documented. This signature is common to other rodents with similar latitudinal range such as Z. hudsonius (King et al. 2006), red-backed voles (Myodes gapperiRunck and Cook 2005), long-tailed voles (Microtus longicaudusConroy and Cook 2000a), deer mice (Peromyscus maniculatusDragoo et al. 2006), woodrats (Neotoma cinereaHornsby and Matocq 2011), chipmunks (Tamias amoenus and T. ruficaudusGood and Sullivan 2001; Demboski and Sullivan 2003; Good et al. 2003, 2008), red squirrels (Tamiasciurus hudsonicusArbogast et al. 2001; Wilson et al. 2005), and flying squirrels (Glaucomys sabrinus and G. volansArbogast 1999). Similarly, soricomorphs, such as Sorex cinereus, S. monticolus, and S. palustris (Demboski and Cook 2001, 2003; Himes and Kenagy 2010), and lagomorphs (e.g., Ochotona princepsGalbreath et al. 2009, 2010) also show this latitudinal signature, suggesting that a common set of processes has influenced diversification across these montane organisms.

Characteristic signatures of population expansion following glacial retreat (Lessa et al. 2003, 2010; Hewitt 2004; Excoffier et al. 2009) include minimal haplotype sorting, low nucleotide diversity, lack of equilibrium between mutation—drift and migration—drift, and starlike phylogeny, which were detected for the Boreal clade of jumping mice. Specific results of tests of neutrality (Table 1) are suggestive of demographic instability since the Last Glacial Maximum. Ancestors of the Boreal clade likely expanded northward in the Holocene due to glacial retreat (Graham et al. 1996; McGill et al. 2005) and in this case, a few closely related haplotypes typify populations ranging from Wyoming (Yellowstone National Park) northward to south-coastal Alaska (Fig. 3), probably the most recently colonized region. At that northern limit, populations separated by 100 km share haplotypes. Post-Pleistocene glacial retreat and the signature of genetic expansion (Table 1) suggest smaller ancestral population size. Ancestral populations were likely restricted to refugia, as documented for other mammals (Waltari et al. 2007; Sommer and Zachos 2009). The newly formed populations at higher latitudes in Canada and southeastern Alaska potentially originated from a source in the south, as reflected by minimal differentiation of haplotypes among populations that span this large area. With the warming climate and retreating glaciers, jumping mice populations likely closely tracked newly available habitats (Hewitt 2004). There is no signal of an isolated refugium in southeastern Alaska as proposed by Jones (1981) for jumping mice and hypothesized for other mammals (Cook et al. 2006).

We documented prolonged isolation (Arenas et al. 2012) for the southem clades, which reflects expected patterns of complete haplotype sorting, deep genetic divergence across the landscape (Fig. 3), high nucleotide and haplotype diversity, and, in general, mutation—drift and migration—drift equilibrium (except Southem Sierra; Table 1). The Southern Sierra clade, however, reflects a significant departure from neutrality, suggestive of a smaller ancestral population (e.g., bottleneck) and a more complex history than our simple hypotheses (ephemeral versus persistent). Genetic footprints of population expansion documented for the Boreal and Southem Sierra clades are likely due to different mechanisms (e.g., latitudinal expansion versus elevational expansion), but this comparison will require more detailed analyses of paleoenvironments and expanded sampling. Elevation fluctuations and concordant genetic signatures have been documented for alpine plants, pika, and woodrats (Matocq 2002b; DeChaine and Martin 2005b; Galbreath et al. 2009, 2010; Beever et al. 2010) across western North America. However, the magnitude of range shifts may differ for species living in montane environments (elevational shifts) versus more homogeneous environments (latitudinal shifts—Guralnick 2007) and these altemative signatures should be explored further (Parmesan 2006; Rubidge et al. 2012). Shifts in elevation and latitude are projected to correspond to changing temperatures (Walther et al. 2002; Parmesan 2006; Petit et al. 2008) and several species in the Sierras have declined over the last century, including jumping mice (Moritz et al. 2008). Other species have shown recent extirpations, such as alpine pika (Galbreath et al. 2009, 2010; Beever et al. 2010), with declines also common elsewhere (DeChaine and Martin 2005a; Knowles and Richards 2005; Albach et al. 2006; Haubrich and Schmitt 2007) including jumping mice (Frey and Malaney 2009).

Niche overlap.—The Sierran jumping mice and Z. trinotatus overlap in niche space more than either does with other lineages of Z. princeps (Fig. 5). This overlap may mirror their close evolutionary relationship; however, overlap may simply reflect the spatial proximity and ecological similarity of the 2 regions. Still, there is significantly more niche divergence between lineages than expected by chance based on pseudoreplicates of background niche space. Because organisms can shift niche preferences through time (Hadly et al. 2009; Peterson 2011), the roles of niche conservatism (Wiens 2004; Wiens and Graham 2005; Warren et al. 2008) or niche divergence (Raxworthy et al. 2007; Rissler and Apodaca 2007) in speciation warrants further exploration. Both have been shown to operate at various temporal and spatial scales in Mexican jays (AphelocomaMcCormack et al. 2010), common kingsnakes (Lampropeltis getulaPryon and Burbrink 2009), and deer mice (P. maniculatusKalkvik et al. 2011). Jumping mice are presumed to have diverged in allopatry based on coalescent simulations (see below; Fig. 4) with presumably both niche conservatism and divergent selection playing roles in the evolution of western jumping mice.

Taxonomic implications.—Westem jumping mice represent a more complex taxonomic assemblage than previously documented. Deep molecular divergence discovered among southern populations of the western jumping mouse, including paraphyly with respect to Z. trinotatus, alters our understanding of species limits in this group (Krutzsch 1954; Hall 1981; Holden and Musser 2005). Hall (1981) recognized 11 subspecies of Z. princeps and 4 subspecies of Z. trinotatus following Krutzsch (1954) review of morphological characters. Our molecular and niche assessment suggests that the initial alignment of the Sierra Nevada populations (Elliot 1898; Preble 1899; Howell 1920) close to Z. trinotatus is appropriate. Gene trees based on mtDNA may not always reflect species limits due to historical mitochondrial introgression (Good et al. 2008; Runck et al. 2009), but in this case, the independent nuclear perspectives corroborate mtDNA and demonstrate the need for revision of species limits in western jumping mice. A comprehensive reevaluation of morphological variation across nominal Z. princeps and Z. trinotatus coupled with development of additional nuclear markers and exploration of finer scale niche variation might provide clarity on the spatiotemporal aspects of diversification.

Historical-biogeographic patterns.—Phylogeographic stracture in the Sierran jumping mice appears to reflect longterm sustained faunal isolation, north-south division of lineages, and elevational shifts with warming climates (Moritz et al. 2008). Other vertebrates in the Sierra Nevada also show a pronounced north-south split, such as wood rats (Matocq 2002a, 2002b; Matocq and Murphy 2007; Matocq et al. 2007), salamanders (Wake 1997), and newts (Tan and Wake 1995). A concordant signature among several species may reflect the influence of glaciers and pluvial lakes formed during the Pleistocene that impeded gene flow (James et al. 2002; Gillespie and Zehfuss 2004). Further evaluation of shifts in elevation through glacial cycles, in combination with comparative assessments of temporal and spatial congruence in lineage diversification across codistributed taxa, is needed for the region.

In general, southem clades of Z. princeps demonstrate strong phylogeographic structure that reflects long periods of isolation without mixing of lineages. There were no haplotypes shared among geographic regions (Fig. 3) with molecular signatures (Table 1) indicative of long-term segregation during the Last Glacial Maximum in a series of isolated refugiai areas across western North America. Coalescent simulations reject an admixture hypothesis but not the persistent allopatric hypothesis (Fig. 4). Jumping mice lineages exhibit higher levels of mtDNA divergence than documented for many other sisterspecies comparisons in mammals (Baker and Bradley 2006). Further, segregation may have persisted over multiple glacial cycles. Multilocus data, coupled with fossil calibration and relaxed molecular clocks (Drummond et al. 2006; Heled and Drummond 2010), have been used to establish initial isolation events in birds and mammals in western North America (McCormack et al. 2011; Reid et al. 2012). Independent lines of evidence suggest there is a common process of allopatric divergence, with historical vicariance via intervening xeric environments, responsible for phylogeographic signatures in DNA and niches among codistributed taxa (Sullivan et al. 2000; Arbogast and Kenagy 2001; Zink 2002; Carstens et al. 2005). Implications of a common signature suggest shared biogeographic processes (Gutiérrez-García and Vázquez Dominguez 2011; Ronquist and Sanmartin 2011) at the community level. Thus, projected climate change and potential shifts in distribution may have more profound (community level) effects than previously considered (Thomas et al. 2004; Moritz et al. 2008; Ackerly et al. 2010).

In conclusion, deep molecular divergence within Z. princeps is accentuated over the southem portion of its current distribution. The wide latitudinal range of Z. princeps provides future opportunities to test hypotheses of incipient speciation using multilocus models and coalescent techniques (Lessa et al. 2003; Carstens et al. 2005). Refinement of the persistent allopatric hypothesis includes testing among vicariant speciation models but serves as a working hypothesis to explore concerted signatures among codistributed species. Examination of these preliminary data suggests that geographic separation between southern lineages has been a dominant and persistent force shaping divergence within Z. princeps and presumably sympatric mammals. Whether these vicariant signatures are suggestive of a common process that is spatially and temporally shared across codistributed mammals, versus simply idiosyncratic responses to fluctuating climate, should be explored.

Supporting Information

Supporting Information S1.—Specimens examined. Found at DOI: 10.1664/12-MAMM-A-006.S1

Acknowledgments

The following institutions provided specimens: University of Alaska Museum of the North, Fairbanks (UAM); Museum of Southwestern Biology, University of New Mexico, Albuquerque (MSB); and Museum of Vertebrate Zoology, University of California, Berkeley (MVZ). Several analyses were completed with the University of Alaska, Fairbanks Life Science Informatics facility (http://biotech.inbre.alaska.edu/portal/). We thank J. Dragoo, F. Torrez-Perez, and B. Truett for assistance with laboratory work; E. Waltari for Excel macro development to down-sample records; members of the Center for Understanding Evolutionary Relationships of Various Organisms (CUERVO) molecular laboratory and S. O. MacDonald for thoughtful discussions; plus 2 reviewers for critical suggestions on a previous version of this manuscript. The project was partially based on support and specimens from the Grinnell Resurvey Project (Yosemite Fund, National Geographic Society grant 8190-07, and National Science Foundation DEB 0640859), Island Surveys to Learn about Endemic Species (ISLES; United States Department of Agriculture Forest Service), and the Beringian Coevolution Project (National Science Foundation 0415668).

Footnotes

  • Associate Editor was Samantha M. Wisely.

Literature Cited

View Abstract