Integrative species delimitation of the widespread North American jumping mice (Zapodinae)

Molecular Phylogenetics and Evolution 114 (2017) 137–152

Contents lists available at ScienceDirect

Molecular Phylogenetics and Evolution journal homepage: www.elsevier.com/locate/ympev

Integrative species delimitation of the widespread North American jumping mice (Zapodinae) Jason L. Malaney a,b,⇑, John R. Demboski c, Joseph A. Cook b a

Department of Biology, Austin Peay State University, Clarksville, TN 37044, USA Museum of Southwestern Biology, University of New Mexico, Albuquerque, NM 87131, USA c Department of Zoology, Denver Museum of Nature & Science, Denver, CO 80205, USA b

a r t i c l e

i n f o

Article history: Received 9 March 2017 Revised 12 May 2017 Accepted 1 June 2017

Keywords: Bayes factors Cryptic diversity Napaeozapus Species distribution model Species tree Zapus

a b s t r a c t Delimiting species can be challenging, but is a key step for the critical examination of evolutionary history and for prioritizing conservation efforts. Because systematic relationships are often determined iteratively using tests based on taxonomy, such methods can fail to detect cryptic variation and result in biased conclusions. Conversely, discovery-based approaches provide a powerful way to define operational taxonomic units and test species boundaries. We compare both approaches (taxonomy-based delimitation – TBD and discovery-based delimitation – DBD) within North American jumping mice (Zapodinae) using broad sampling, multilocus analyses, and ecological tests. This group diversified through the dynamic glacial-interglacial periods of the Quaternary and phylogeographic tests reveal 28 lineages that correspond poorly with current taxonomy (4 species, 32 nominal subspecies). However, neither the 4-species or 28-lineage hypotheses are optimal for species-level classification. Rather, information theoretic approaches (Bayes Factors) indicate a 15-species hypothesis is best for characterizing genetic variation in this group, with subsequent iterative pairwise ecological tests failing to confirm four species pairs. Taken together, evolutionary and ecological tests capture divergence among 11 putative species that, if upheld by additional tests, will lead to taxonomic revision and reevaluation of conservation plans. Ó 2017 Elsevier Inc. All rights reserved.

1. Introduction Understanding how divergence ultimately leads to speciation is a central question in evolutionary biology. Application of molecular phylogenetic approaches routinely reveals patterns of differentiation that improve our collective understanding of both the origin and timing of divergence, and when coupled with independent data (e.g., ecological variation) can expose underlying diversification processes (Crandall et al., 2000; Wiley and Lieberman, 2011). Refining our understanding of evolutionary relationships and the ecogeographic limits of taxa provides a requisite context for other fields (Bernardo, 2011; Riddle and Hafner, 1999), particularly conservation (Haig et al., 2006). Species delimitation approaches aim to refine our understanding of diversification processes (Camargo and Sites, 2013; Carstens et al., 2013; Jackson et al., 2016). Nevertheless, because most species likely evolve in

⇑ Corresponding author at: P.O. Box 4718, Department of Biology, Austin Peay State University, Clarksville, TN 37044, USA. E-mail addresses: [email protected] (J.L. Malaney), [email protected] (J.R. Demboski), [email protected] (J.A. Cook). http://dx.doi.org/10.1016/j.ympev.2017.06.001 1055-7903/Ó 2017 Elsevier Inc. All rights reserved.

allopatry, a complication for any delimitation technique is whether sufficient evolutionary and ecological changes have accrued to allow species detection (Coyne and Orr, 2004; de Queiroz, 2007; Mayr, 1963). For example, species long isolated by a barrier can accumulate neutrally evolving genetic differences, but because environments may remain similar, few ecological or even morphological differences may be apparent. Conversely, recently diverged species can rapidly adapt to distinctive localized conditions, often resulting in ecological and morphological divergence, but few neutrally evolving genetic changes. Because of these different evolutionary signatures, delimitations based solely on a single class of characters may fail to detect speciation. Therefore, species delimitation techniques that integrate across multiple datasets should be favored (Carstens et al., 2013; Sukumaran and Knowles, 2017). Traditionally taxonomists have faced a choice among a variety of competing conceptual approaches to either delimit or describe species (Camargo and Sites, 2013; Sites and Marshall, 2004; Wiens, 2007). Species delimitation has experienced resurgence due to new integrative approaches (Carstens et al., 2013; Edwards and Knowles, 2014; Hope et al., 2016) that are based on hypothesis tests that better capture spatiotemporal signatures of

138

J.L. Malaney et al. / Molecular Phylogenetics and Evolution 114 (2017) 137–152

evolutionary history (Grummer et al., 2013; Padial et al., 2010; Yang and Rannala, 2010). Analyses that integrate diverse datasets and apply explicit tests are expected to be most powerful when coupled with coalescence methods (Carstens and Dewey, 2010; Fujita et al., 2012). There are several species delimitation approaches, therefore the choice of approach for characterizing operational taxonomic units (OTUs) for testing for the signal of speciation is nontrivial (Carstens et al., 2013; Sukumaran and Knowles, 2017). Approaches to species delimitation include two broad techniques. Unguided approaches attempt to jointly estimate OTUs and test phylogenetic relationships, but they remain problematic and have received less attention despite active development (Ence and Carstens, 2011; Jones et al., 2015; O’Meara, 2009; Yang and Rannala, 2014). Conversely, guided approaches require OTUs from the outset and are now common (Edwards and Knowles, 2014; Toprak et al., 2016). Unguided approaches will likely supplant guided approaches, but there currently remains no consensus over an optimal technique. Critically, guided approaches that depend on existing taxonomy from the outset require careful consideration, because a flawed taxonomic hypothesis can introduce undesirable bias (Olave et al., 2014). Guided approaches use two broadly defined methods to establish OTUs (Ence and Carstens, 2011). The first uses taxonomy usually drawn from an existing premise of morphologically identified species and subspecies. Taxonomy-based delimitation (TBD), also called validation-based approaches (O’Meara, 2009), assesses taxonomic units by applying validation tests of described variation. The second method uses discovery-based delimitation (DBD), also called substructure detection methods (Rannala, 2015). When applying DBD, a set of OTUs or ‘candidate species’ is identified using genetic clustering or assignment tests, thereby providing a starting point to compare and evaluate alternative hypotheses of systematic relationships (Camargo and Sites, 2013; Edwards and Knowles, 2014; Leache and Fujita, 2010). In some cases, taxonomic and discovery-based approaches produce similar results, but in situations where taxonomy is suspect, DBD may provide alternative perspectives on species limits. Both TBD and DBD require robust tests using multiple datasets rather than relying on a singular dataset (Fujita et al., 2012; Rannala, 2015). Morphology-based taxonomy often forms a priori hypotheses in species delimitations in vertebrates. Mammalian systems have been exemplar models for species-tree approaches and existing taxonomy frequently guides sampling (Ence and Carstens, 2010; Heled and Drummond, 2010). While taxonomy-based tests may represent a natural progression from earlier studies (i.e., validation of described morphological variation), they may also violate assumptions of some analyses, especially when paraphyly or gene flow are present. Furthermore, tests predicated on previously described variation often miss cryptic variation (Leache et al., 2009; Pyron and Burbrink, 2009). Tests predicated solely on morphology may best reflect phenotypic responses to environments, but may also fail to provide the optimal basis for exploring evolutionary history or establishing taxonomy (Gotthard and Nylin, 1995; Mayr, 1956; Simpson, 1951), although there are exceptions (Hoekstra et al., 2005; Patton and Smith, 1994). Consequently, delimitation approaches that rely exclusively on described phenotypic variation (i.e., TBD) should be validated with independent data and analyses (Carstens et al., 2013; Fujita et al., 2012; Rannala, 2015). In this study, we unite signatures of spatial, genetic, and ecological divergence to compare and evaluate alternative hypotheses of species limits in a continentally-distributed clade of rodents. More specifically, we apply guided approaches to species-delimitation and evaluate whether taxonomy optimally reflects evolutionary history. We generate hypotheses that originate from phylogeographic signatures (DBD) that we compare to morphologically-

based taxonomy (TBD). After contrasting DBD and TBD approaches using genetic data, we then apply iterative ecological tests to assess genetic conclusions. Taken together, we present a new taxonomic hypothesis that better reflects phylogeny, biogeography, ecology, and has immediate implications for conservation.

1.1. Study system Jumping mice (Dipodidae, Zapodinae) are broadly distributed (Fig. 1) from northern Alaska and Canada, south to the montane areas of the American Southwest and Appalachians (Hall, 1981; Krutzsch, 1954). Zapodines include two genera, four species, and 32 subspecies that are spread over 30 North American ecoregions (SF1, SF2; (Holden and Musser, 2005; Krutzsch, 1954). Such a diverse set of environments provides an excellent system to assess how physical barriers, geographic isolation, and environmental variability interact to influence divergence (Pigot et al., 2010). Phylogeographic studies have uncovered cryptic variation and paraphyly between the western jumping mouse (Zapus princeps) and Pacific jumping mouse (Z. trinotatus (Malaney et al., 2013). Multiple Pacific coastal forest and montane-associated lineages have significantly deeper divergences than observed within the more broadly distributed grassland-associated meadow jumping mouse (Z. hudsonius; (King et al., 2006; Malaney and Cook, 2013; Malaney et al., 2012; Ramey et al., 2005). Similarly, the woodland jumping mouse (Napaeozapus insignis) reflects a deep allopatric divergence associated with the temperate forests of eastern North America (Malaney and Cook, 2013). When combined, at least four, major phylogeographic patterns exist across the wide distribution of jumping mice (Malaney and Cook, 2013) that are generally shared with several other widespread mammals, suggesting a common response to climate fluctuations in northern-temperate mammals (Malaney, 2012). A consensus among jumping mice phylogeographic studies is that there is deep and often cryptic variation within species, and evidence that evolutionary lineages are discordant with taxonomy (Himes and Kenagy, 2013; King et al., 2006; Malaney et al., 2013; Malaney and Cook, 2013; Malaney et al., 2012; Ramey et al., 2005). Consequently, testing taxonomic hypotheses against phylogeographic hypotheses is crucial for optimally characterizing variation and defining conservation priorities. More specifically, this clade is ideal for testing DBD and TBD approaches because these phylogeographic signals (Himes and Kenagy, 2013; Malaney et al., 2013; Malaney and Cook, 2013) challenge existing classifications (Hall, 1981; Holden and Musser, 2005; Krutzsch, 1954).

2. Methods 2.1. Generalized approach We used molecular sequences and ecological variables to establish and test candidate species (discovery-based delimitation, DBD) and to test described variation (taxonomy-based delimitation, TBD), using broad geographic sampling that includes all recognized species and subspecies. For the discovery-based hypotheses, we apply four approaches expected to capture the spectrum of evolutionary and ecological variation and that reasonably sample patterns of speciation (Sukumaran and Knowles, 2017). We then evaluate hypotheses generated from both TBD and DBD using multilocus species-trees and Bayes Factor tests to identify the hypothesis that best reflects evolutionarily significant diversity (Grummer et al., 2013). Finally, we iteratively apply pairwise ecological tests as an approximation of niche conservatism/divergence among candidate species (sister taxa) predicted by the optimal hypothesis.

J.L. Malaney et al. / Molecular Phylogenetics and Evolution 114 (2017) 137–152

139

Fig. 1. Taxonomy-based delimitation (TBD) for North American Zapodines. Geographic distributions (left) and species-tree phylogeny (right) of the four species revised by Krutzsch (1954; K54:4), but that accommodates paraphyly among Z. princeps and Z. trinotatus. Symbols reflect sampling localities of 117 samples applied in all species-tree analyses.

For TBD, we applied OTUs from either recognized (Hall, 1981; Krutzsch, 1954) or proposed (Jones, 1981) taxonomies using species and subspecies designations. Establishing OTUs using DBD is a more involved process. First, we used a network-based approach based on mtDNA because phylogenetic networks have served as the basis for establishing conservation units and validating some subspecies within this system (Crandall and Marshall, 2007; King et al., 2006; Ramey et al., 2005, 2006; Vignieri et al., 2006). Network-based approaches can be problematic for a variety of reasons, so, second, we identified significant mitochondrial clade structure using an evolutionary model (Fujisawa and Barraclough, 2013; Reid and Carstens, 2012). Third, because mtDNA-based inference is considered limited (Ballard and Whitlock, 2004; Hickerson et al., 2006; Knowles and Carstens, 2007), nuDNA variation is often favored in species delimitations, so we applied genetic clustering approach (Satler et al., 2013). Finally, delimitations based on a single class of characters (i.e., genetic) can be biased, so we used an approach that jointly considers spatial, genetic, and ecological variation to identify a set of candidate species (Edwards and Knowles, 2014). Statistical species delimitation using multilocus DNA coupled with coalescence model-selection is a rapidly evolving field (Edwards and Knowles, 2014; Grummer et al., 2013; Olave et al., 2014). In general, species delimitation requires three steps: (1) identify OTUs, (2) determine the phylogeny relating putative groups, and (3) compare alternative hypotheses. We used DBD and TBD to complete the first step, applied species-tree approaches for the second step, and used a pair of model-selection techniques during the final step. For model selection, we applied Bayes Factors (2lnBf) to discriminate between alternative hypotheses. Both simulation tests and empirical data indicate these approaches are robust, even in systems with introgression or stochastic lineage sorting (Grummer et al., 2013). 2.2. Taxon sampling We genetically sampled populations of all recognized genera (Hall, 1981; Krutzsch, 1954), and proposed (Jones, 1981) species

and subspecies of North American jumping mice (genera: Napaeozapus and Zapus), obtaining 431 samples from both natural history collections (N = 379) and fieldwork (N = 52) that targeted type localities of species and subspecies. We added 488 GenBank records representing phylogeographic and cryptic genetic variation within subspecies (Himes and Kenagy, 2013; King et al., 2006; Ramey et al., 2005). All samples were identified to species, with subspecies designations added using a combination of morphology and geographic distributions (Hall, 1981; Himes and Kenagy, 2013; Malaney et al., 2013; Malaney and Cook, 2013). Sample localities and GenBank numbers are provided in Appendix 1, with more detailed specimen information for most samples on Arctos (http://arctos.database.museum/). 2.3. Genetic datasets To characterize genetic variation across jumping mice, we used multilocus approaches by including sequences from five genes, including 441 mitochondrial cytochrome b sequences (cytb; 1140 bp) from our previous work (Malaney et al., 2013, 2012; Malaney and Cook, 2013), and another 322 sequences (1006 bp) from King et al. (2006). We used nuDNA from a subset (N = 117) of representative individuals including Apolipoprotein B (APOB; 347 bp), Breast Cancer Susceptibility (BCRA1; 824 bp), Glucocerebrosidase (GBA; 346 bp), and the myosin heavy chain polypeptide (MYH6; 267 bp) gene all generated from a previous study (Malaney et al., 2013; Malaney and Cook, 2013). We identified polymorphic sites using the IUPAC nucleic acid code. We applied an unphased framework because phased datasets have high computational costs and require excessive parameterization, but may not improve phylogenetic estimates among allopatrically diverged organisms when using species-tree approaches (Kubatko et al., 2011; Kubatko and Gibbs, 2010). For analyses, we partitioned the genetic data in three ways. First, we analyzed the 117 samples with both mtDNA and nuDNA data (five genes) to explore alternate signals derived from each locus. Second, we applied coalescent-based species-tree analyses to those 117 samples, but that partitioned terminal taxa differently for hypothesis tests. Because mtDNA

140

J.L. Malaney et al. / Molecular Phylogenetics and Evolution 114 (2017) 137–152

Because cluster-based approaches may not reflect evolutionary relationships and tree-based methods may fail to reveal reticulate evolution in recent divergences (Posada and Crandall, 2001), we used phylogenetic networks to identify discrete evolutionary units for the mtDNA dataset. We conducted network analyses with tcs v. 1.21 (Clement et al., 2000). This analysis implements the statistical-parsimony approach of Templeton et al. (1992) and may offer a better exploration of incipient lineages and intraspecific diversity (Avise, 2008). Finally, we applied Gaussian clustering that jointly considers spatial, genetic, and ecological variation to identify a set of candidate species and has demonstrated high performance in species delimitations (Edwards and Knowles, 2014; Hausdorf and Hennig, 2010). We generated individual pairwise Gower distances (Gower, 1971) for spatial and environmental data using the daisy, cluster, and vegan packages in R (Maechler et al., 2016; Oksanen et al., 2008; R Core Team, 2014). For the spatial distance, we calculated Euclidian distance between genetically sampled locations (latitude and longitude). We calculated environmental distance using 20 environmental variables from the WorldClim database (www.worldclim.org), and included the Ecoregion (Level II) obtained from the Environmental Protection Agency (www.epa.gov). Multivariate analyses (McGarigal et al., 2000) can be sensitive to correlation so we removed excessively correlated environmental variables using Pearson’s correlation coefficient. We calculated individual pairwise genetic distances using MEGA v. 5.0 (Tamura et al., 2011) and applied corresponding evolutionary model corrections (Table 1) to account for multiple substitutions (Felsenstein, 2001). To account for differences in substitution rates, we partitioned individual loci by mean pairwise distances and individual distances were averaged across loci. We standardized distance matrices for each dataset (spatial, genetic, environmental) using NMDS with isoMDS using the MASS package (Venables and Ripley, 2002). We followed Hausdorf and Hennig (2010) by retaining four dimensions for each dataset. Stress values for NMDS were generally <10%, with absolute maximum of 20% and we concatenated the four NMDS dimensions for each dataset (spatial + genetic + ecological).

can bias species-tree analyses, we conducted a second set of species-tree analyses by excluding mtDNA from all 117 samples. Finally, we used the complete mtDNA dataset (N = 763) to identify lineage-specific distributions during niche-based tests because that dataset represents the densest geographic sampling (Fig. 1). 2.4. Genetic analyses To reconstruct evolutionary relationships, we conducted phylogenetic analyses of individual genes using Bayesian inference. First, we subjected each aligned dataset to alternative models of sequence evolution in jModelTest (Posada, 2008), and applied Bayesian Information Criterion (Posada and Buckley, 2004) to determine the best-fit nucleotide substitution model (Table 1). For each gene, we conducted phylogenetic reconstructions with MrBayes v.3.2.1 (Huelsenbeck and Ronquist, 2005; Ronquist et al., 2011) that were initiated with random trees and ran with four chains (default heating values) for five million generations, sampling every 5 k, with duplicate analyses to ensure convergence. We assessed the target distribution of likelihood scores by evaluating split frequencies of the standard deviation and examined optimal parameter estimates in Tracer (Drummond and Rambaut, 2007). We also assessed convergence by heuristically inspecting the diagnostic plots of the Metropolis-coupled Markov chain Monte Carlo (MC3). We identified nodal support (posterior probability – PP) using the consensus of the residual trees with the first 10% of iterations discarded as burn-in (Huelsenbeck and Imennov, 2002) and depicted final trees with FigTree v. 1.3.1 (SF3; (Rambaut and Drummond, 2009). We used an implementation of the general mixed Yulecoalescent model (GMYC) to account for uncertainty in tree topology and branch lengths in our single-locus mtDNA dataset. This method may best reflect the spectrum of diversification (Pons et al., 2006), has been widely used in species delimitations (e.g., (Esselstyn et al., 2012; Papadopoulou et al., 2008; Satler et al., 2013), and is robust to simplifying assumptions, including the principal factor to delimitations, the mean population size relative to divergence time (Fujisawa and Barraclough, 2013). We used the R package bGMYC (Reid and Carstens, 2012) to sample across the posterior distribution of gene trees using 100 evenly sampled trees estimated from MrBayes that we rescaled to be ultrametric using the ape package (Paradis et al., 2004). In general, default bGMYC settings were retained, except the py2 parameter that we altered to 1.2 (program documentation) and the initial number of species was set to half the number of tips. We conducted analyses for 20 k generations with burn-in set to 10k and sampling every 200 generations. Because nuclear genomes are chimeric objects composed of fragments of potentially differing origins and histories, we estimated genetic clusters of nuDNA genetic variation that applied Bayesian clustering (BC), which may capture complex patterns of subdivision better than other metrics (e.g., Fst) (Pritchard et al., 2000). We applied BC using the program Structure v. 2.3.4 using default settings and genotypes derived from nuDNA sequences (Pritchard et al., 2000, 2010).

2.4.1. Species tree reconstruction At the species level, single-gene analyses often indicate a lack of monophyly and topologies may vary among loci (Carstens and Knowles, 2007; Knowles and Kubatko, 2010; Maddison, 1997), so we estimated species-tree phylogenies using a multilocus approach. At a minimum, we sampled at least one individual representing each OTU for all loci but our sampling generally included sequences from at least three individuals per OTU across all genes (Camargo et al., 2012). We conducted species-tree reconstruction using the ⁄BEAST approach in BEAST v.1.8.3 (Drummond and Rambaut, 2007) using the CIPREs web service. When applying fossil calibrations, Malaney et al. (2012) found that strict and relaxed clock methods did not significantly change inferences in this system. In this study, fossils are unsuitable for most hypotheses complicating comparisons. Therefore, because our primary interest is relative diversification rates and

Table 1 Summary of nucleotide variation for multi-locus analyses including the abbreviated name of each locus, number of sites (bp), exon-intron spanning regions, number of individuals (N), mean number of nucleotide differences (k), number of segregating sites (S), total mutations (g), number of haplotypes (h), number of parsimony informative sites (PI), and the optimal model of DNA evolution. Locus

bp

Region

N

k

S

g

h

Gap

PI

Model

cytb APOB BRCA1 GBA MYH7

1140 367 789 346 267

Entire ex26 ex11 ex4-5 ex36-37

117 117 117 117 117

155.21 2.05 13.52 3.41 4.01

451 23 180 24 29

607 29 216 26 37

101 35 50 22 49

– 3 4 3 2

426 19 167 24 26

GTR+I+C HKY HKY+C HKY JC+I

141

J.L. Malaney et al. / Molecular Phylogenetics and Evolution 114 (2017) 137–152

comparisons across hypotheses (Table 2), we enforced a strict clock for the mtDNA (rate 0.05) for analyses and allowed estimated substitution rates (i.e., relaxed local molecular clocks) for nuDNA (Drummond and Suchard, 2010; Drummond et al., 2006; Welch and Bromham, 2005). For each species-tree analysis, we conducted a minimum of 50 m generations sampled every 50 k with convergence diagnostics examined in Tracer v. 1.5. We identified the single topology that best represents the posterior distribution of all plausible trees using TreeAnnotator. Summary parameters included burn-in of 10% and 0.5 posterior probability limits, with mean node heights for divergence estimates. We depicted the final maximum clade credibility species trees using FigTree. Because mtDNA often has strong signal that may bias speciestree estimates (Knowles, 2010; Knowles and Kubatko, 2010), we also reconstructed each hypothesis (TBD & DBD) by excluding mtDNA. For these analyses, we enforced a strict clock for just the BCRA gene (rate 0.01, based on estimates from Malaney and Cook, 2013), allowed for relaxed local molecular clocks for the remaining loci, and retained all other parameters to maintain consistency with the mtDNA + nuDNA analyses.

2.4.2. Systematic hypotheses We compared seven alternative hypotheses that utilize the multilocus dataset (N = 117) for species-tree analyses and Bayes Factors model selection (Table 2). For the first three hypotheses, we applied TBD. The first hypothesis (K54:4) is predicated on existing taxonomy (N. insignis, Z. hudsonius, Z. princeps, Z. trinotatus) last revised by Krutzsch (1954), and summarized in Holden and Musser (2005). Because some subspecies are frequently considered independent evolutionary units and have been the focus of conservation (e.g., Z. h. preblei and Z. h. luteus), we considered the 32 subspecies as OTUs in the second hypothesis (K54:32). For the third taxonomy-based hypothesis (J81:2), we considered only two species as suggested by Jones (1981; N. insignis and uniting all Zapus). We contrasted the three TBD hypotheses with four DBD hypotheses (Table 2). Two hypotheses were based on solely on mtDNA. First, we assessed whether 28 mtDNA lineages optimally represent evolutionarily significant variation using a Genealogical approach. Next, we implemented a Coalescent approach (Fujisawa and Barraclough, 2013; Zhang et al., 2011, 2013) that has been widely used in species delimitations (Esselstyn et al., 2012; Papadopoulou et al., 2008; Satler et al., 2013). Independent nuDNA loci are often favored for DNA-based species delimitations (Carstens et al., 2013; Leache and Fujita, 2010; Olave et al., 2014), so we applied a Chimeric approach to generate an estimate of the optimal nuclear genetic clusters of jumping mice using only nuDNA. Finally, because singular character types (i.e., genetic) can fail, we applied an Integrative (spatial + ecological + genetic)

approach using Gaussian clustering (Edwards and Knowles, 2014). Following species-tree analyses for each hypothesis, we conducted model comparisons to identify the hypothesis that best explained the data considering marginal likelihood estimates of competing hypotheses using both Path and Stepping Stone sampling (Grummer et al., 2013). We ranked scores for each hypothesis with respect to the hypothesis with the lowest marginal likelihood. We then considered the strength of support for the leading hypothesis compared to competing hypotheses (Table 2). We accepted the optimal hypothesis when both model comparisons were >10 suggesting strong support for the leading topology in pairwise comparisons. In the case where these criteria were <10, we selected the model with the fewest parameters using the likelihood score penalized by the number of parameters in each model (Burnham and Anderson, 2002, 2004). 2.5. Species distribution modeling We used SDMs to identify important niche variables and compare ecological variation among candidate species within a phylogenetic context (McCormack et al., 2010) and then test for nichebased variation and ecological divergence (Araújo and Peterson, 2012). We applied niche-based tests using a two-step approach. First, we reconstructed contemporary distributions by uniting jumping mice occurrence records from natural history museums with several environmental variables presumed to be important in their ecological divergence. Second, we applied pairwise comparisons of putative taxa to the niche identity test to assess whether they occupy an identical distribution of environmental variables (Warren et al., 2008, 2010). 2.5.1. Variables and extent Temperature and precipitation variables are important predictors of jumping mice distributions (Malaney and Cook, 2013; Malaney et al., 2012), and are presumably important in the Grinnellian niche space (Soberón, 2007). Therefore, we used 19 abiotic bioclimatic variables from across the study region (Hijmans et al., 2005), plus elevation from the WorldClim database (http:// www.worldclim.org) each at 2.5 min resolution or 4 km2. The extent of the study area can strongly influence predictive models of species distributions (Anderson and Raza, 2010; Barve et al., 2011; Graham and Hijmans, 2006). Therefore, we restricted the variables to a region that represents the known distribution of each species and added a 250 km buffer to define the study extent. Models can be influenced by highly correlated variables (Rissler and Apodaca, 2007), but there remains no consensus of the optimal approach to incorporate correlation into species distribution models. Nevertheless, we randomly distributed 10,000 points across

Table 2 Model comparison among seven alternative hypotheses of species limits for North American jumping mice (Zapodinae). Hypotheses were derived from either discovery-based delimitation (DBD) or taxonomy-based delimitation (TBD). Two Bayes Factor tests were used to evaluate log marginal likelihoods and compare and evaluate hypotheses using model selection including stepping-stone sampling (SS) and path sampling (PS) using MCMC in BEAST. The hypothesis with the optimal score is indicated in bold, however for the nuDNA three models were close (D < 10.0) so we chose the model with the fewest parameters (Integrative: 15 candidate species) that was then used for ecological validation tests. * In both the K54 hypotheses, OTUs that accommodated paraphyly (Malaney et al., 2013) were used rather than enforcing a strict taxonomy prior. All Loci -LN

nuDNA SS

DSS

DPS

DSS

DPS

Hypothesis

Description

DBD

Genealogical Coalescent Integrative Chimeric

28 species - network 21 species - GMYC 15 species -Gaussian 5 species - Structure

16920.2 16915.8 16901.7 16932.0

17543.4 17568.9 17592.5 17915.8

0 25.58 49.15 372.52

17546.9 17566.4 17590.9 17912.9

0 19.52 44.01 366.02

7528.7 7531.4 7525.7 7540.3

7937.7 7938.6 7942.2 8150.5

0 0.89 4.49 212.74

7935.5 7936.9 7940.4 8148.7

0 1.41 4.92 213.23

TBD

K54:4 K54:32 J81:2

4 species* - Krutzch 1954 32 subspecies* - Krutzch 1954 2 species – Jones, 1981

16930.6 16955.1 16926.9

18035.8 17696.8 18380.3

492.42 153.44 836.98

18030.8 17692.0 18372.5

483.88 145.11 825.55

7551.7 7590.2 7550.9

8234.1 7971.4 8489.5

296.38 33.68 551.78

8231.8 7968.8 8485.2

296.37 33.30 549.76

PS

-LN

SS

PS

142

J.L. Malaney et al. / Molecular Phylogenetics and Evolution 114 (2017) 137–152

study areas using ArcGIS and extracted values for each variable and conducted a Pearson’s product-moment correlation coefficient (r) using the rcorr routine in the hmisc package of R (R Core Team, 2014). We identified and removed variables with r |>0.80| that can complicate interpretations of niche models (Rissler and Apodaca, 2007). 2.5.2. Occurrences For SDMs, we downloaded localities for each species (19,328 voucher-based records of all jumping mice) from VertNet.org (April 2016) and georeferenced 4000 of these specimens using Biogeomancer. We then screened records for potentially erroneous data and removed those with >5 km radius error. We re-coded species and subspecies assignments based on recognized taxonomy and assigned proposed candidate species designations. Presence-only modeling can be prone to sampling biases that may result in a model representing an over-fit distribution (Kramer-Schadt et al., 2013; Reddy and Davalos, 2003; Veloz, 2009), so we used a thinning approach that discards redundant sampling locations within a 10 km radius yielding 1072 down-sampled records used for modeling. 2.5.3. Correlative models We used environmental variables and occurrence records to reconstruct SDMs by applying Maximum Entropy (Elith et al., 2011; Graham et al., 2008; Merow et al., 2013). These approaches can offset potential biases due to implementing presence-only correlative modeling. All models were constructed using the program Maxent version 3.3.3a (Elith et al., 2006; Phillips et al., 2006) with 20 replicate runs and random background sampling within study extents (Phillips et al., 2009). Species-specific parameter tuning often enhances model performance (Anderson and Gonzalez, 2011), so we conducted a preliminary analysis of model selection (Warren and Seifert, 2011) to assess the optimization of the b parameter using ENMTools v.1.3 (Warren et al., 2010). During modeling, we used localities with genetic information as training datasets to provide validation for the predicted distribution in the final models whenever possible. When insufficient numbers of genetic samples were available, we randomly reserved 20% of localities for training. To explore anomalous model behavior that may be indicative of locality biases or inappropriate background sampling, we examined standard deviation across all replicates. We assessed model performance using the area under the operating characteristic curve (AUC) which is a threshold-independent measure ranging from 0 to 1 (Peterson et al., 2008) and useful for evaluating models constructed from presence/background samples (Fielding and Bell, 1997; Hanley and McNeil, 1982). Maxent produces a set of continuous surfaces and we used the median of replicated models to represent relative climatic suitability (Richmond et al., 2010). There are a variety of approaches for establishing thresholds (Jimènez-Valverde and Lobo, 2007; Liu et al., 2005), and we converted the continuous surface to a binary distribution using a 90% threshold criterion (Pearson et al., 2007) that represented suitable and unsuitable regions. Final models were assessed by comparing the distributions of genetically verified individuals to reconstructed distributions (Fig. 4). 2.5.4. Tests of ecological divergence Closely related taxa tend to share niche space but frequently fail to occupy identical ecological conditions, thus ecological divergence often represents a continuum ranging from niche divergence to conservatism (Warren et al., 2008; Wiens and Graham, 2005). Pairwise comparisons are frequently used to assess ecological divergence among closely related species. For hypothesis tests of ecological divergence we applied niche identity tests using Schoener’s D (Schoener, 1968; Warren et al., 2008). We used niche

identity tests to assess if candidate species occupy an identical distribution of ecological conditions (one-tailed test), which is based on randomization of sample points and reconstructing an expected distribution with values ranging between 0.0 (discordant) and 1.0 (identical). We used ENMTools to assess niche identity and conduct randomization tests using 100 replicates. 3. Results 3.1. Genetic data & phylogenetic patterns We obtained sequences from five genes with complete taxon sampling representing all 32 subspecies and denser infraspecific geographic sampling for the mtDNA dataset, except subspecies with narrow ranges in degraded habitats (e.g., Z. p. curtatus and Z. p. chrysogenys; see Malaney, 2012), where only a single locality was sampled. Sequence variation and genetic diversity measures between nominal species and subspecies was relatively high (Malaney and Cook, 2013). We applied the datasets in two ways. First, we used 1006 bp of mtDNA (N = 763) to define phylogeographic and genealogical variants and confirm niche-based tests. For reconstructing evolutionary relationships, we used 1140 bp of mtDNA and 1784 bp of nuDNA for multilocus species-tree analyses (N = 117). Optimal DNA substitution models indicated different models of evolution for each locus and that loci were unlinked (Table 1). Reconstructions of evolutionary history of mtDNA using a GMYC revealed 21 phylogenetically distinct clades with the 95% confidence interval (CI) capturing between 18 and 25: two within N. insignis, six in Z. hudsonius, seven in Z. princeps (one cryptic), and three in Z. trinotatus. Network analyses of mtDNA generally corresponded with the phylogenetic trees but detected a finer level of variation with seven additional divisions for a total of 28 lineages detected from 289 haplotypes (Appendix 1; Fig. 2). Inspection of mtDNA trees reveals strong support for 25% (8/32) of the morphologically based subspecies. However, these trees also reveal paraphyletic relationships between most subspecies and even between species (Malaney et al., 2013). Gene trees for the nuDNA reflect varying levels of resolution (SF4) that generally correspond to deeper mtDNA clades. We also identify previously undetected cryptic variation within Z. princeps and poor support for most (21/32) recognized subspecies. None of the nuDNA trees reflected strong support for individual subspecies. 3.2. Species-tree reconstruction & systematic hypotheses Phylogenetic reconstructions of individual genes are often incongruent following recent divergence. We detected incongruence among species and candidate species for individual nuDNA trees and we also detected paraphyletic relationships among nominal jumping mice species (SF4). With certain accommodations, species-tree analyses can account for incongruence (Knowles and Carstens, 2007; Knowles and Kubatko, 2010; Maddison, 1997) and when we applied accommodations to species tree analyses for known paraphyly among nominal Z. princeps and Z. trinotatus we detect improved performance by ⁄BEAST (Table 2). Regardless of gene-tree species-tree incongruence, we examined strong posterior distributions of all parameters of interest in species-tree analyses, except those that failed to account for paraphyly (SF2). In all species-tree analyses, we detected strong support for the initial late-Pliocene divergence of the group with subsequent diversification occurring throughout the Pleistocene (Figs. 2, 4, 5). In general, divergence times coincide with glacial-interglacial transitions, but 95% highest posterior densities (HPD) for some nodes transcend glacial/interglacial events (Fig. 5). In all species-tree

J.L. Malaney et al. / Molecular Phylogenetics and Evolution 114 (2017) 137–152

143

Fig. 2. Alternative species delimitation hypotheses for North American jumping mice. The top row represents the recognized taxonomy (Carleton and Musser, 2005; Hall, 1981) revised by Krutzsch (1954) that includes four species (K54:4) and 32 subspecies (K54:32). The second row represents an alternative taxonomy proposed by Jones (1981; J81:2). This study presents four novel discovery-based hypotheses. At the finest scale, 28 phylogeographic lineages (Genealogical) were used to generate the speciestree phylogeny derived from five genes (phylogeny). Stars above lineages represent conservation units, including lineage 11 (USFWS Endangered: Z. h. preblei plus Z. h. tenellus and Z. h. alascensis) and lineage 12 (USFWS Endandered: Z. h. luteus). Note the cryptic diversity for lineage 14 and paraphyletic relationships for lineages 23 and 24 (Z. p. pacificus; Malaney et al., 2013). Within the phylogeny, open circles represent high nodal support (posterior probabilities 0.95) and deeper nodes with poor support (nodes X and Y) correspond to periods of rapid diversification or deep coalescent events (node Z). 2 M and 4 M indicate 2 and 4 million years ago respectively. Inset photo: USFWS Endangered, New Mexico Meadow Jumping Mouse (Z. h. luteus; Lineage 12) from Bosque del Apache National Wildlife Refuge, NM, USA. Photo credit: Greg Wright.

reconstructions, except J81:2, we detected three deep nodes (Figs. 1, 2, 4, 5). This includes an initial divergence 4.5 MYA between the genera Napaeozapus and Zapus (Fig. 2) and then divergence at roughly 3 MYA and at 2 MYA, resulting in four strongly supported clades (Clades I–IV). However, at least 14 strongly supported nodes were detected within the four clades that may be indicative of speciation. Taken together, we detected, with strong

support (95% PP), at least 14 divergences that have occurred within the last 1 Million years. In contrast, using the genealogical 28-lineage phylogeny (Fig. 2), we failed to detect strong support (<95% PP) for nine nodes, typically near tips or during periods of rapid divergence. Among the deeper nodes (those 500 kya or older), ⁄BEAST analyses identified poor support for three nodes associated with short branches

144

J.L. Malaney et al. / Molecular Phylogenetics and Evolution 114 (2017) 137–152

Fig. 3. AADDCCBBNiche-based geographic distributions of 15 candidate species of North American jumping mice. The four panels correspond to the four deep clades (current species-level lineages, see Figs. 2 and 4). In all distribution models, binary thresholds were applied (see methods) and colors correspond with proposed taxonomy (Fig. 4) including Clade I: A = Red, B = Navy; Clade II: C-E = Green, F = Brown; Clade III: G = Charcoal, H&I = Yellow, J = Orange, K = Blue, Clade IV: L&M = Spruce, N = Grey, O = Magenta. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

associated with rapid divergence and large effective population sizes (Figs. 2 and 4). This included node X within nominal Z. hudsonius, an ancestor of eight mtDNA lineages (4–11, Fig. 2) that diverged 500 kya and is now associated with grasslands and northern latitudes. Second, node Y within nominal Z. princeps, an ancestor that diverged 1.50 MYA and gave rise to populations that now occupy the Okanogan Valley and northern Rocky Mountains (14–20, Fig. 3). Finally, node Z within nominal Z. trinotatus represents the ancestor of five lineages (23–27, Fig. 3) that diverged 750 kya and descendant populations are now found in the Coastal and Sierra Nevada ranges. Alternative species-tree phylogenies (TBD + DBD) but detected differences among hypotheses indicative of the ability of the models to account for observed genetic and coalescent variation. For example, while similar deep divergences were identified, species trees often differed in estimates of coalescent events and effective population sizes. When we excluded mtDNA from species tree reconstructions, performance was indistinguishable (DSS  10 and DPS  10) between the 15-candidate species, the 21-species, and the 28-species models (Table 2), so we retained the model with fewest parameters (15-species) for correlative modeling and iterative ecological tests.

3.3. Correlative models Correlative modeling performed well when interpolating climatic suitability for jumping mice with all mean AUC values >0.80 and low variation across replicates (ST1). Generally, correlative models predicted geographically restricted distributions that were a product of alternative sets of variables and, in general, the importance of environmental variables was different for predicting the distributions of each candidate species. We detected limited range overlap among candidate species, especially among montane associated candidate species. Exceptions include the overlapping distributions (Fig. 3) of candidate species of jumping mice in the Sierra Nevada range (L and M) and several candidate species (C, D, E) occupying the Great Plains to the high latitudes. In contrast, candidate species associated with barriers typically reflect high range overlap that may be a result of geographic isolation. For example, candidate species on either side of the St. Lawrence River (lineages A and B) and Columbia River (lineages N and O) had high predicted-range overlap. Hypothesis tests of ecological divergence using niche identity detected significant ecological differences between most candidate species (Fig. 4, SF5). However, niche identity tests detected a lack

J.L. Malaney et al. / Molecular Phylogenetics and Evolution 114 (2017) 137–152

145

Fig. 4. Pairwise tests of ecological interchangeability among candidate species of North American jumping mice. Candidate species that failed the validation tests include C, D, and E; H from I; and L from M resulting in 11 putative species proposed for taxonomic revision. Lower Panel: Species-tree phylogeny of the corresponding 15 candidate species with branch widths scaled to estimated Ne (dmv median) and open circles representing high nodal support (0.95 posterior probability).

of ecological differentiation among four pairs of candidate species including C vs E, D vs E, H vs I, and L vs M.

candidate species revealed four pairs with insignificant niche differences, yielding 11 putative species that optimally capture both evolutionary and ecological variation.

4. Discussion 4.1. Alternative systematic hypotheses Characterizing the signal of speciation can be difficult. In this study, we contrasted seven hypotheses of systematic relationships within the North American jumping mice. Taxonomic hypotheses included one proposed and two existing taxonomies (subspecies and species) along with four novel discovery-based hypotheses. For all seven hypotheses, we treated terminal taxa as putative species and applied multilocus species-tree analyses to a broad geographic sampling that ensured complete taxon sampling (Hall, 1981; Holden and Musser, 2005; Krutzsch, 1954). We evaluated support using Bayes Factor tests that initially detected a 28species hypothesis as superior, but information theoretic techniques revealed a 15-species hypothesis as the optimal representation of genetic variation. Iterative ecological tests among the 15

Formal evaluation of alternative hypotheses of divergence is a key step for understanding the mechanisms responsible for generating and maintaining diversity. Species-tree analyses simultaneously accommodate gene-tree and species-tree uncertainty under a common set of assumptions (Heled and Drummond, 2010) such as the presence of phylogenetic discord (Knowles and Kubatko, 2010). In general, we found deep divergences concordant across mtDNA and nuDNA with less concordance, resulting in lower nodal posterior probabilities, near the tips or during rapid diversification (Fig. 2, SF4). Likelihood-based tree calibration methods such a ⁄BEAST require robust OTUs from the outset. Because existing jumping

146

J.L. Malaney et al. / Molecular Phylogenetics and Evolution 114 (2017) 137–152

Fig. 5. Proposed taxonomy and evolutionary relationships for North American jumping mice, which includes 11 species (colored boxes) that began diversifying during the Pliocene (4.6 MYA). Below are major geologic epics with shading corresponding to dominant glacial-interglacial periods (light and dark respectively) throughout the Pliocene (blue), Pleistocene (grey), and Holocene (black). The background topologram depicts all possible hierarchies of individual gene trees encountered during the specietree analysis superimposed with the maximum clade credibility species tree (10% burnin) with bars (purple) signifying the 95% credible intervals around the mean node ages and numbers with corresponding nodal posterior probabilities. 2 M and 4 M indicate 2 and 4 million ago years respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

mice taxonomy may be suspect, we applied an approach to compare hypotheses predicated on taxonomy versus novel discovery-based OTUs. Since Bayes Factor tests revealed that discovery-based approaches generally outperformed taxonomybased approaches (Table 2), the current taxonomy, largely based on skull and skeletal characters (Jones, 1981; Krutzsch, 1954), appears to improperly represent diversity. Consequently, taxonomic hypotheses warrant further discussion because for at least six decades (Krutzsch, 1954) they have guided validation tests (e.g., King et al., 2006; Ramey et al., 2005) and are the basis of conservation (USFWS, 2010, 2014) with significant economic implications (Brook et al., 2003; Crifasi, 2007; Industrial Economics, 2010; Johnson, 2004). Poor performance of the K54:4 hypothesis is potentially the product of two interrelated issues, paraphyly and cryptic variation. Paraphyly between Z. princeps and Z. trinotatus (Himes and Kenagy, 2013; Malaney et al., 2013) could be the result of stochastic coalescent sorting, gene flow, or improperly defined species limits. Malaney et al. (2013) used coalescent simulations coupled with niche and range overlap tests to show that jumping mice occupying the Sierra Nevada (currently assigned to Z. princeps) are more closely related to Z. trinotatus. We assigned individuals to terminal taxa using both strict taxonomic assignment (ignoring paraphyly) and by accommodating paraphyly (Fig. 2). Nevertheless, even by accommodating paraphyly, the K54:4 model failed to yield significant improvement over alternative hypotheses, which is likely the product of coalescent events that remain unaccounted for within each putative species. Thus, alternative hypotheses that can accommodate this cryptic variation displayed improved performance (e.g., Genealogical, Table 2). For the K54:32 hypothesis, we also assigned individual samples to subspecies by both accommodating and ignoring paraphyly (SF2). We expected this approach to be problematic because of poor match between genetic variation and morphologically

described variation for some subspecies. While in a few cases (i.e., Z. h. luteus, Z. p. chrysogenys) we detect concordance between subspecies descriptions and genetic/ecological variation, in others (e.g., Z. p. utahensis, N. i. abietorum), subspecies are products of multiple evolutionary lineages (Fig. 2, Table 3). We also detect instances of mtDNA haplotype sharing among subspecies, often indicative of a recent history. For example, the widespread Northern lineage is composed of the few genetic differences between Z. h. alascensis, Z. h. preblei, and Z. h. tenellus (Malaney and Cook, 2013). Furthermore, some subspecies are composed of individuals found to be part of two OTUs. For example, both Z. h. americanus and Z. h. canadensis have individuals of both candidate species C & D (Table 3). Finally, in at least one case (i.e., Okanogan Valley and eastern slopes of the Cascade Range), cryptic variation remains undescribed (i.e., lineage 14, candidate species G). Taken together, the K54:32 model performed poorly compared to most DBD hypotheses, but is perhaps a reasonable approximation of morphological differences representing ecogeographic variation. Nevertheless, the K54:32 hypothesis presents an incomplete prediction of evolutionary divergence, a common theme among mammals (Ashton et al., 2000; Mace, 2004; Ryder, 1986). We found that DBD approaches that defined alternative hypotheses of OTUs from the outset were nearly always superior to TBD hypotheses. Two DBD hypotheses represent either marginal improvement (Coalescent) or underperform (Chimeric) when compared to other DBD hypotheses, so we limit our discussion to the Genealogical and Integrative hypotheses. We consider the implications of the Genealogical hypothesis due to its superior performance (Table 2), and because mtDNA-based delimitations are often favored in mammalian taxonomy (Baker and Bradley, 2006) plus analogous approaches were applied to define conservation units in this system (King et al., 2006; Ramey et al., 2005). The Genealogical hypothesis performed exceptionally well in all Bayes Factor tests (Table 2). However, single-locus delimitations

J.L. Malaney et al. / Molecular Phylogenetics and Evolution 114 (2017) 137–152 Table 3 Summary of current and proposed taxonomies including four species and 32 subspecies, plus 15 candidate species that performed best in Bayes Factor assessments and used in ecological validation tests. Results suggest eleven species proposed for taxonomic revision. Nomenclature for 10 of the species are resurrected from previously described variants (Hall, 1981; Krutzsch, 1954; Preble, 1899) and one is proposed as a new species (sp nov). *Zapus princeps pacificus is paraphyletic and more closely related to Z. trinotatus (Figs. 2 and 4, Malaney et al., 2013), but presented here as described by Krutzsch (1954). Current taxonomy Species

Subspecies

Candidate

Proposed taxonomy

Napaeozapus insignis

abietorum part frutectanus saquenayensis abietorum part insignis roanensis

A

Napaeozapus abietorum (Preble, 1899)

acadicus

C

americanus canadensis ladas alascensis campestris hudsonius intermedius preblei tenellus luteus pallidus

C&E C&E C D D D D D D F F

NONE cinereus curtatus oregonus utahensis part chysogenys princeps saltator utahensis part idahoensis kootenayensis minor pacificus*

G H H H&I H

Zapus okanoganensis sp nov Zapus oregonus (Preble, 1899)

J J K K

Zapus princeps (Allen, 1893)

K K K L&M

eureka

N

montanus orarius trinotatus part trinotatus part

N N N

Zapus hudsonius

Zapus princeps

Zapus trinotatus

A A B

Napaeozapus insignis (Miller, 1893)

B B

O

Zapus hudsonius (Zimmermann, 1780)

Zapus luteus (Miller, 1911)

Zapus saltator (Allen, 1893)

Zapus pacificus (Elliot, 1898; Merriam, 1897) Zapus orarius (Merriam, 1897)

147

Exceptional performance of the Genealogical hypothesis is likely a product of coalescent signal captured by the hypothesized terminal taxa, even among the nuDNA, that individually fail to accurately capture all mtDNA lineages (SF3). Incongruence between mtDNA and nuDNA is either the product of incomplete lineage sorting or gene flow. Characterizing sources of incongruence was not the focus of the current study, but there is growing evidence, especially among mammals of divergence in the face of gene flow (Morales et al., 2016; Nosil, 2008; Sullivan et al., 2014). Because the 28-species hypothesis has >50% of coalescent events occurring during the last 100,000 years (Fig. 2) with 82% (23/28) of the lineages seemingly allopatric, incomplete lineage sorting may be the dominant source of incongruence. Consequently, there has likely been insufficient time for nuDNA differences to accumulate. More detailed geographic and genomic sampling is needed to more rigorously assess areas of sympatry among lineages and decouple the source of incongruence. However, despite superior performance of the Genealogical hypothesis (28-species model) for both mtDNA + nuDNA and nuDNA species trees, we did not recognize this model as optimal because it presents only marginal improvement over two other discoverybased hypotheses (Coalescent & Integrative; Table 2) when applying only nuDNA and information theoretic Bayes Factor tests. The Integrative hypothesis (15-species) was inferior to the Genealogical hypothesis (28-species) when mtDNA was included but performed similarly to other hypotheses with mtDNA excluded (Table 2). Consequently, we recognize the less parameterized Integrative hypothesis (15-species) as the optimal characterization of speciation among jumping mice. While gene tree topologies are generally congruent the mtDNA and individual nuDNA coalescent signals differ somewhat across genes, a result common in species-tree analyses (Kubatko et al., 2011; Leache et al., 2009). Therefore, the signal detected by mtDNA tests may be capturing lineage specific signal incompletely captured by nuDNA and perhaps reflects geographic (sub)structure but not necessarily speciation. Protracted speciation model simulations suggest that delimitations frequently detect structure, not speciation (Sukumaran and Knowles, 2017). For jumping mice, some lineages may experience episodic divergence with periods of secondary contact allowing admixture that perhaps results in patterns expected by the protracted model (e.g., across Z. hudsonius). Additional tests are needed, but our results may suggest that Bayes Factors accurately detect signals of speciation under a protracted speciation model. 4.2. TBD vs DBD

Zapus trinotatus (Rhoads, 1894)

can be problematic (Carstens et al., 2013; Edwards and Knowles, 2014; Olave et al., 2014). Therefore, we estimated species trees by using approaches that both include and exclude mtDNA to explore both genomes. During all analyses we allowed gene trees to embed within species trees, thereby tolerating gene tree–species tree discordance and accommodating lack of monophyly for all gene trees (Knowles and Carstens, 2007). Lack of monophyly among gene trees is often interpreted as incomplete lineage sorting, but such patterns still contain information about species relationships (Knowles, 2010; Kubatko et al., 2011). In this system, we detect incomplete sorting among all gene trees (SF3). Nevertheless, the Genealogical hypothesis resulted in optimal performance whether we included mtDNA in species-tree analyses (Fig. 2) or treated the nuDNA as a validation test of the mtDNA generated hypothesis (SF4, Table 2).

Our results highlight two broad conclusions related to taxonomy. First, while TBD is commonly applied in validation of taxonomic assignments, it can result in flawed tests. Therefore, second, comparisons between TBD and DBD are valuable for assessing alternative hypotheses of systematic relationships, especially in polytypic systems, and should likely be used more frequently. TBD is often favored for validating variation but can be subjective if criteria are poorly defined (Carstens et al., 2013). Jumping mice taxonomy based on four species with 32 subspecies present an unsuitable prior for species delimitations or even iterative validation tests. In contrast, DBD approaches avoid issues that arise by applying taxonomy a priori in species delimitations. Integrative approaches are frequently favored because singular sets of characters (i.e., genetic) can yield an incomplete perspective of divergence. In jumping mice, the optimal hypothesis identified by model selection was derived by applying an Integrative approach (Table 2). Some studies recognize species using the finest set of divisions such as monophyly of mtDNA (Camargo et al., 2012;

148

J.L. Malaney et al. / Molecular Phylogenetics and Evolution 114 (2017) 137–152

Esselstyn et al., 2012; Kubatko et al., 2011; Leache and Fujita, 2010). We detect strong support using the Genealogical approach that detected 28 evolutionarily significant lineages. Nevertheless, approaches that fail to delimit the finest variation are likely optimal for stabilizing taxonomy (Rosell et al., 2010; Satler et al., 2013; Sukumaran and Knowles, 2017). Although some may advocate recognizing 28 species of jumping mice based on significant genetic differences among the mtDNA dataset, Bayes Factors failed to detect significant differences between the Genealogical, Coalescent, and Integrative approaches when using nuDNA (Table 2). Furthermore, there are strong ecological similarities among some closely related lineages. When combined, organellar differences recently accumulated do not necessarily correspond to nuclear differences or ecological variation. Alternative partitioning schemes of terminal taxa may influence delimitations, especially among putative species with low posterior support or in regions with overlapping ranges or in secondary contact, and these should be explored in more detail. Discrepancies between DBD and TBD approaches may reflect alternate axes of variation. For example, morphology can be influenced by ecophenotypic variation (i.e., morphological response to environment) that may not closely mirror evolutionary relationships, as is often the case in rodents (Caumul and Polly, 2005). Alternatively, morphological variation may have been poorly characterized during initial descriptions when several subspecies were characterized using few specimens. For example, initial description of Z. h. preblei was based on four individuals (Krutzsch, 1954). Our data are unsuitable to distinguish whether ecophenotypic variation influenced taxonomy or if morphological variation was poorly characterized. Because we detect at least 15 species using genetic approaches, additional tests are needed to assess if geographic variation in morphological characters best match a 15-species model or the taxonomy of Krutzsch (1954). If morphology best matches taxonomic models, then North American jumping mice may be morphologically constrained despite significant genetic and ecological differences across their range. 4.3. Ecological variation We integrated ecological variation among jumping mice into analyses in two ways. First, we included environmental information to detect candidate species using Gaussian clustering that united signals of spatial, genetic, and environmental variation. Second, we reconstructed SDMs and applied niche-based tests to the 15species hypothesis (Fig. 4). Taken together, candidate species of jumping mice have significant ecological differences, which is perhaps not surprising considering they occupy at least 30 level III ecoregions (SF1-2). However, perhaps most surprising is the significant ecological differences detected between pairs of candidate species not previously incorporated into taxonomy. While additional multivariate tests are likely needed, in cases of apparent lack of ecological interchangeability detected by niche equivalency tests, we infer that divergence stems from a variety of processes (Pyron et al., 2015; Shafer and Wolf, 2013). Because both environmental and physical heterogeneity is high in this system, we expect ecophysical allopatry (i.e., environmental and geographic isolation) is likely responsible for observed ecological divergence, whereas tests that failed to reject the hypothesis of ecological interchangeability may reflect niche conservatism (Pyron et al., 2015; Wiens et al., 2010; Wiens and Graham, 2005). Reconstructed distributions of jumping mice are reflective of ecological differences among candidate species and we detect a general lack of spatial overlap between the predicted distributions and low ecological interchangeability between candidate species using pairwise tests (Fig. 3). Additional phylogeographic sampling is needed for all candidate species to validate and codify these

hypothesized distributions. We detected seven cases of low interchangeability and four cases where ecological differences were insignificant (Fig. 4, SF5), suggesting a variety of factors may contribute to jumping mice diversification including conservatism/constraint for the latter four cases (Pyron et al., 2015). Candidate species L and M both occupy the Sierra Nevada may be separated by a physiological barrier (Malaney et al., 2013). Both candidate species are predicted by a similar set of environmental variables (ST1), with precipitation seasonality (Bio15) having the greatest contribution, permutation importance, the most information (highest gain), and the greatest information not present in other variables (decreased gain). Relative contributions of other variables differed slightly between these candidate species (ST1). Consequently, seasonal variation in precipitation is especially important for predicting the distribution of both candidate species. Because year-to-year precipitation can be extremely variable in the region, with expected increases in amplitude (Mann and Gleick, 2015), additional tests should evaluate the potential impacts of climate change on local populations. Ecological divergence was not supported between candidate species C and D, C and E, and D and E. When combined, candidate species C + D + E represents the widest ranging species encompassing a bicoastal distribution (>100 degrees of longitude) with genetic signals of post-glacial demographic expansion (Malaney and Cook, 2013). Similar to tests in this study, Malaney and Cook (2013) detected few genetic or ecological differences between populations occupying the Front Range of Colorado and Wyoming and those found further north (e.g., Alaska). Limited genetic or ecological differences over wide geographic areas indicates that additional phylogeographic and phylogenetic tests are needed to clarify evolutionary history (e.g., Clade II, Figs. 2–4) and decouple factors contributing to the lack of ecological differences across this wide range. Environmental variables useful for predicting the distributions of these three candidate species (C, D, E) are similar, with Bio1 (mean annual temperature), Bio9 (temperature of driest quarter), Bio12 (annual precipitation), and Bio15 (precipitation seasonality) dominating models (ST1). We detected few ecological differences between candidate species H and I from the Intermountain West (Fig. 3), which may reflect the strong imprint of Holocene climatic change driving upslope contraction of previously connected grasslands and forests across this large geographic area. Candidate species H occurs in Uinta Mountains whereas candidate species I occupies several mountain rages within and adjacent to the Great Basin. These populations diverged during the last interglacial and now occupy a variety of environmental conditions across this arid region. This high variability of environmental conditions may have contributed to our inability to detect significant environmental differences predicted by genetic differences of these candidate species. Across the Intermountain West, higher elevations are often mesic and forested with interspersed valleys typically xeric and lacking forests. Background sampling during modeling included both mountains and valleys potentially affecting our ability to distinguish among ecological conditions. Because of few ecological differences between these candidate species, despite significant signals of genetic differences, we fail to recognize them as distinct species until further tests are conducted. However, in situations where we detect both significant genetic and ecological differences (Fig. 4), we infer speciation has occurred suggesting a revised taxonomy should include 11 species of jumping mice. 4.4. Taxonomic recommendations The perspective of 11 putative species of jumping mice (Fig. 5) is predicated on (1) a reasonable hypothesis generated using an integrative discovery-based technique, (2) tested using multilocus

J.L. Malaney et al. / Molecular Phylogenetics and Evolution 114 (2017) 137–152

species-tree analyses, (3) compared to reasonable alternative hypotheses using Bayes Factors and information theoretic approaches, and that were (4) assessed using pairwise ecological tests. Taken together, this new perspective better captures the deep evolutionary history (seven of the ten divergences are at least 1 million years old) and ecological differences observed across jumping mice. Future work should assess range boundaries, characterize biogeographic history, and test species limits predicated on assignment of individuals and levels of potential gene flow. Comparisons with morphological variation are needed before formal taxonomic revision, but we tentatively recommend a revised nomenclature (Fig. 5, Table 3). 4.5. Sampling effects A century of taxonomic and systematic studies of jumping mice highlights how alternative species concepts and disparate delimitation criteria can produce discordant classifications, which raise two important issues related to sampling. First, previous efforts (Ramey et al., 2005; King et al., 2006) shared a common approach that relied solely on an existing taxonomic framework to guide sampling and test hypotheses. In this study, we identified that use of an existing taxonomy can be problematic. Systematic treatments reasonably begin with the assumption that described morphological variation accurately predicts genetic variation, with species and subspecies designations often used as OTUs for validation tests. Previous validation studies applied subspecies designations to assess genetic variation among populations of Z. hudsonius occupying the Front Range of Colorado and Wyoming, and part of the Great Plains (King et al., 2006; Ramey et al., 2005). Those studies arrived at opposite conclusions despite using similar data, prompting extensive debate over taxonomic treatments (Crandall and Marshall, 2007; Cronin, 2007; Vignieri et al., 2006), statistical inference (Brosi and Biber, 2009), and the role of politics in science (Carolan, 2008; Crifasi, 2007; Lackey, 2007; Scott et al., 2007). We sampled all subspecies, often from type localities, but as shown by our analyses, at least some subspecies designations of jumping mice are tenuous, so we recommend subspecies be carefully re-evaluated before they form the basis of research or management actions. Second, previous studies of jumping mice failed to include complete taxon sampling and assumed adjacent subspecies were closely related and so missed important genetic variation and failed to account for biogeographic history (Malaney and Cook, 2013). Biogeographic history, in particular, offers an important foundation for conservation (Moritz and Potter, 2013). For example, populations from the Front Range of Colorado and Wyoming were the source for populations now found in British Columbia and Alaska, due to the rapid, northward expansion of jumping mice during the Holocene, a common pattern in high latitude mammals (Lessa et al., 2003; Malaney, 2012). Therefore, populations along the Front Range have a biogeographic history largely independent from adjacent subspecies occupying the southern Great Plains, a finding that compromises conclusions of previous work, but highlights that dense geographic sampling is crucial for taxonomic clarity. 4.6. Implications for conservation Conservation efforts concentrated solely on a few geographically restricted populations may miss others of greater conservation value that may be experiencing elevated threat (Malaney and Cook, 2013). The expanded evolutionary and ecological perspective presented in this study should help to restructure conservation priorities for jumping mice. In general, some federally listed populations are minimally divergent, while others with restricted

149

distributions and deep divergence and are perhaps imperiled without formal protection. Ideally, conservation efforts should capture evolutionary diversity, however formalized protections (USFWS, 2010, 2014) for jumping mice are focused within a single clade (clade II, Fig. 2) that diverged during the late Pleistocene. The two federally listed subspecies (i.e., Z. h. luteus and Z. h. preblei) represent <5% of observed variation. In contrast, clades III and IV are distributed across the western North America (candidate species G-O, Figs. 2 and 3), have experienced episodic diversification throughout the Pleistocene, and now account for a majority of the genetic (87%) and ecological variation observed within Zapus. Furthermore, some of these candidate species have isolated populations that may be in decline. The current taxonomic perspective of North American jumping mice predates a rigorous understanding of phylogenetic relationships, demographic processes, and evolutionary history yet all are important considerations for conservation (Gutiérrez and Helgen, 2013). We have demonstrated that an emphasis on morphological variation can fail to adequately capture evolutionary history and ecological differences, and ultimately miss significant evolutionary units that may be critical to conservation efforts. Moreover, this study reinforces the value of applying approaches that more completely capture the signal of evolutionary diversification, and importantly, yields an improved path forward for characterizing and conserving biodiversity. Acknowledgements We recognize the support and feedback from multiple individuals including the CUERVO, Matocq, and Feldman labs. We are particularly grateful to Kayce Bell, Bryan McLean, Jack Sullivan, Enrique Lessa, and two anonymous reviewers for providing critical feedback on earlier versions of this manuscript. We are especially appreciative for the access to voucher specimens and loaned materials from various museums and acknowledge these essential resources offered by specimen-based collections and various field biologists whose work is too frequently underappreciated. Without those resources and efforts, this study would not have been possible. Funding for this research was provided by the USFWS and NSF (DEB 1258010). Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.ympev.2017.06. 001. References Allen, J.A., 1893. List of Mammals collected by Mr. Charles P. Rowley in the San Juan Region of Colorado, New Mexico, and Utah, with Descriptions of New Species. Bull. Amer. Mus. v, pp. 69–84. Anderson, R.P., Gonzalez, I.J., 2011. Species-specific tuning increases robustness to sampling bias in model of species distributions: An implementation with Maxent. Ecol. Model. 222, 2796–2811. Anderson, R.P., Raza, A., 2010. The effect of the extent of the study region on GIS models of species geographic distributions and estimates of niche evolution: preliminary tests with montane rodents (genus Nephelomys) in Venezuela. J. Biogeogr. 37, 1378–1393. Araújo, M.B., Peterson, A.T., 2012. Uses and misuses of bioclimatic envelope modeling. Ecology 93, 1527–1539. Ashton, K.G., Tracy, M.C., de Queiroz, A., 2000. Is Bergmann’s rule valid for mammals? Am. Nat. 156, 390–415. Avise, J.C., 2008. Three ambitious (and rather unorthodox) assignments for the field of biodiversity genetics. Proc. Natl. Acad. Sci. U.S.A. 105, 11564–11570. Baker, R.J., Bradley, R.D., 2006. Speciation in mammals and the genetic species concept. J. Mammal. 87, 643–662. Ballard, J.W.O., Whitlock, M.C., 2004. The incomplete natural history of mitochondria. Mol. Ecol. 13, 729–744. Barve, N., Barve, V., Jimenez-Valverde, A., Lira-Noriega, A., Maher, S.P., Peterson, A.T., Soberón, J., Villalobos, F., 2011. The crucial role of the accessible area in

150

J.L. Malaney et al. / Molecular Phylogenetics and Evolution 114 (2017) 137–152

ecological niche modeling and species distribution modeling. Ecol. Model. 222, 1810–1819. Bernardo, J., 2011. A Critical Appraisal of the Meaning and Diagnosability of Cryptic Evolutionary Diversity, and Its Implications for Conservation in the Face of Climate Change. Climate Change, Ecology and Systematics. Cambridge University Press, Cambridge. Brook, A., Zint, M., De Young, R., 2003. Landowners’ responses to an endangered species act listing and implications for encouraging conservation. Conserv. Biol. 17, 1638–1649. Brosi, B.J., Biber, E.G., 2009. Statistical inference, Type II error, and decision making under the US Endangered Species Act. Front. Ecol. Environ. 7, 487–494. Burnham, K.P., Anderson, D.R., 2002. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. Springer. Burnham, K.P., Anderson, D.R., 2004. Multimodel inference understanding AIC and BIC in model selection. Sociol. Methods Res. 33, 261–304. Camargo, A., Avila, L.J., Morando, M., Sites, J.W.J., 2012. Accuracy and precision of species trees: effects of locus, individual, and base pair sampling on inference of species trees in lizards of the Liolaemus darwinii group (Squamata, Liolaemidae). Syst. Biol. 61, 272–288. Camargo, A., Sites, J.W., 2013. Species delimitation: a decade after the renaissance. In: Pavlinov, I. (Ed.), The Species Problem – Ongoing Issues. InTech. Carleton, M.D., Musser, G.G., 2005. Order rodentia. In: Wilson, D.E., Reeder, D.M. (Eds.), Mammal Species of the World. The Johns Hopkins University Press, Baltimore, Maryland. Carolan, M.S., 2008. The politics in environmental science: The Endangered Species Act and the Preble’s mouse controversy. Environ. Polit. 17, 449–465. Carstens, B., Dewey, T.A., 2010. Species delimitation using a combined coalescent and information-theoretic approach: an example from North American Myotis bats. Syst. Biol. 59, 400–414. Carstens, B., Pelletier, T.A., Reid, N.M., Satler, J.D., 2013. How to fail at species delimitation. Mol. Ecol. 22, 4369–4383. Carstens, B.C., Knowles, L.L., 2007. Estimating species phylogeny from gene-tree probabilities despite incomplete lineage sorting: An example from Melanoplus grasshoppers. Syst. Biol. 56, 400–411. Caumul, R., Polly, P.D., 2005. Phylogenetic and environmental components of morphological variation: skull, mandible, and molar shape in marmots (Marmota, Rodentia). Evolution 59, 2460–2472. Clement, M., Posada, D., Crandall, K.A., 2000. TCS: a computer program to estimate gene genealogies. Mol. Ecol. 9, 1657–1659. Coyne, J.A., Orr, H.A., 2004. Speciation. Sinauer Associates, Sunderland, Massachusetts. Crandall, K.A., Bininda-Emonds, O.R., Mace, G.M., Wayne, R.K., 2000. Considering evolutionary processes in conservation biology. Trends Ecol. Evol. 15, 290–295. Crandall, K.A., Marshall, J.C., 2007. An assessment of the threatened subspecific status of the Preble’s meadow jumping mouse (Zapus hudsonius preblei) based on current molecular data sets. Genoma LLC: Report For the State of Wyoming, Office of the Attorney General, Woodland Hills, UT, pp. 1–35. Crifasi, R.R., 2007. A subspecies no more? A mouse, its unstable taxonomy, and western riparian resource conflict. Cult. Geograph. 14, 511–535. Cronin, M.A., 2007. The Preble’s meadow jumping mouse: subjective subspecies, advocacy and management. Anim. Conserv. 10, 159–161. de Queiroz, K., 2007. Species concepts and species delimitation. Syst. Biol. 56, 879– 886. Drummond, A., Suchard, M., 2010. Bayesian random local clocks, or one rate to rule them all. BMC Biol. 8, 114. Drummond, A.J., Ho, S.Y., Phillips, M.J., Rambaut, A., 2006. Relaxed phylogenetics and dating with confidence. PLoS Biol. 4, e88. Drummond, A.J., Rambaut, A., 2007. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol. Biol. 7. Edwards, D.L., Knowles, L.L., 2014. Species detection and individual assignment in species delimitation: can integrative data increase efficacy? Proc. Roy. Soc. London B: Biol. Sci. 281, 20132765. Elith, J., Graham, C.H., Anderson, R.P., Dudik, M., Ferrier, S., Guisan, A., Hijmans, R.J., Huettmann, F., Leathwick, J.R., Lehmann, A., Li, J., Lohmann, L.G., Loiselle, B.A., Manion, G., Moritz, C., Nakamura, M., Nakazawa, Y., Overton, J.M., Peterson, A.T., Phillips, S.J., Richardson, K., Scachetti-Pereira, R., Schapire, R.E., Soberón, J., Williams, S., Wisz, M.S., Zimmermann, N.E., 2006. Novel methods improve prediction of species’ distributions from occurrence data. Ecography 29, 129– 151. Elith, J., Phillips, S.J., Hastie, T., Dudík, M., Chee, Y.E., Yates, C.J., 2011. A statistical explanation of MaxEnt for ecologists. Divers. Distrib. 17, 43–57. Elliot, D.G., 1898. Lists of species of mammals principally rodents obtained by W. W. Price, Dr. S. E. Meek, G. R. Cherrie, and E. S. Thompson in the states of Iowa, Wyoming, Montana, Idaho, Nevada, and California with descriptions of new species. Field Columbian Museum Publ. 27, 193–221. Ence, D., Carstens, B., 2010. SpedeSTEM: a rapid and accurate method for species delimitation. Mol. Ecol. Resour. 24, 473–480. Ence, D.D., Carstens, B.C., 2011. SpedeSTEM: a rapid and accurate method for species delimitation. Mol. Ecol. Resour. 11, 473–480. Esselstyn, J.A., Evans, B.J., Sedlock, J.L., Khan, F.A.A., Heaney, L.R., 2012. Single-locus species delimitaiton: a test of the mixed Yule-coalescent model, with an empirical application to Philippine round-leaf bats. Proc. Roy. Soc.ety B-Biol. Sci. 279, 3678–3686. Felsenstein, J., 2001. Taking variation of evolutionary rates between sites into account in inferring phylogenies. J. Mol. Evol. 53, 447–455.

Fielding, A.H., Bell, J.F., 1997. A review of methods for the assessment of prediction errors in conservation presence/absence models. Environ. Conserv. 24, 38–49. Fujisawa, T., Barraclough, T.G., 2013. Delimiting species using single-locus data and the Generalized Mixed Yule Coalescent (GMYC) approach: a revised method and evaluation on simulated datasets. Syst. Biol. 62, 707–724. Fujita, M.K., Leache, A.D., Burbrink, F.T., McGuire, J.A., Moritz, C., 2012. Coalescentbased species delimitaiton in an integrative taxonomy. Trends Ecol. Evol. 27, 480–488. Gotthard, K., Nylin, S., 1995. Adaptive plasticity and plasticity as an adaptation: a selective review of plasticity in animal morphology and life history. Oikos, 3– 17. Gower, J.C., 1971. A general coefficient of similarity and some of its properties. Biometrics, 857–871. Graham, C.H., Elith, J., Hijmans, R.J., Guisan, A., Peterson, A.T., Loiselle, B.A., 2008. The influence of spatial errors in species occurrence data used in distribution models. J. Appl. Ecol. 45, 239–247. Graham, C.H., Hijmans, R.J., 2006. A comparision of methods for mapping species ranges and species richness. Glob. Ecol. Biogeogr. 15, 578–587. Grummer, J.A., Bryson, R.W., Reeder, T.W., 2013. Species delimitation using Bayes factors: simulations and application to the Sceloporus scalaris species group (Squamata: Phrynosomatidae). Syst. Biol., syt069 Gutiérrez, E.E., Helgen, K.M., 2013. Mammalogy: Outdated taxonomy blocks conservation. Nature 495, 314. Haig, S.M., Beever, E.A., Chambers, S.M., Draheim, H.M., Dugger, B.D., Dunham, S., Elliott-Smith, E., Fontaine, J.B., Kesler, D.C., Knaus, B.J., Lopes, I.F., Loschl, P., Mullins, T.D., Sheffield, L.M., 2006. Taxonomic considerations in listing subspecies under the US Endangered Species Act. Conserv. Biol. 20, 1584–1594. Hall, E.R., 1981. The Mammals of North America. John Wiley & Sons, New York, New York. Hanley, J.A., McNeil, B.J., 1982. The meaning and use of the area under a Receiver Operating Characteristic (ROC) curve. Radiology 143, 29–36. Hausdorf, B., Hennig, C., 2010. Species delimitation using dominant and codominant multilocus markers. Syst. Biol., syq039 Heled, J., Drummond, A.J., 2010. Bayesian inference of species trees from multilocus data. Mol. Biol. Evol. 27, 570–580. Hickerson, M.J., Meyer, C.P., Moritz, C., 2006. DNA barcoding will often fail to discover new animal species over broad parameter space. Syst. Biol. 55, 729– 739. Hijmans, R.J., Cameron, S.E., Parra, J.L., Jones, P.G., Jarvis, A., 2005. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978. Himes, C.M., Kenagy, G.J., 2013. Delimiting geographic distributin and populaiton history of jumping mice (Zapus trinotatus and Zapus princeps) in the Pacific Northwest. Northwestern Naturalist 94, 22–34. Hoekstra, H.E., Krenz, J.G., Nachman, M.W., 2005. Local adaptation in the rock pocket mouse (Chaetodipus intermedius): natural selection and phylogenetic history of populations. Heredity 94, 217–228. Holden, M.E., Musser, G.G., 2005. Family dipodidae. In: Wilson, D.E., Reeder, D.M. (Eds.), Mammal Species of the World. The Johns Hopkins University Press, Baltimore, Maryland, pp. 871–893. Hope, A.G., Malaney, J.L., Bell, K.C., Salazar-Miralles, F., Chavez, A.S., Barber, B.R., Cook, J.A., 2016. Revision of widespread red squirrels (genus: Tamiasciurus) highlights the complexity of speciation within North American forests. Mol. Phylogenet. Evol. 100, 170–182. Huelsenbeck, J.P., Imennov, N.S., 2002. Geographic origin of human mitochondrial DNA: Accommodating phylogenetic uncertainty and model comparison. Syst. Biol. 51, 155–165. Huelsenbeck, J.P., Ronquist, F., 2005. Bayesian analysis of molecular evolution using MrBayes. Statist. Methods Mol. Evol., 183–232 Industrial Economics, I., 2010. Economic Analysis of Critical Habitat Designation for Preble’s Meadow Jumping Mouse in Colorad. Final Economic Analysis I, Cambridge, MA. Jackson, N.D., Carstens, B.C., Morales, A.E., O’Meara, B.C., 2016. Species delimitation with gene flow. Syst. Biol., syw117 Jimènez-Valverde, A., Lobo, J.M., 2007. Threshold criteria for conversion of probability of species presence to either-or presence-absence. Acta Oecol. 31, 361–369. Johnson, K., 2004. Debate Swirls Around the Status of a Protected Mouse. The New York Times, New York. Jones, G., Aydin, Z., Oxelman, B., 2015. DISSECT: an assignment-free Bayesian discovery method for species delimitation under the multispecies coalescent. Bioinformatics 31, 991–998. Jones, G.S., 1981. The Systematics and Biology of the Genus Zapus (Mammalia, Rodentia, Zapodidae). The School of Graduate Studies. Indiana State University, Terre Haute, Indiana. King, T.L., Switzer, J.F., Morrison, C.L., Eackles, M.S., Young, C.C., Lubinski, B.A., Cryan, P., 2006. Comprehensive genetic analyses reveal evolutionary distinction of a mouse (Zapus hudsonius preblei) proposed for delisting from the US Endangered Species Act. Mol. Ecol. 15, 4331–4359. Knowles, L.L., 2010. Chapter 10 Sampling strategies for species tree estimation. In: Knowles, L.L., Kubatko, L. (Eds.), Estimating Species Trees: Practical and Theoretical Aspects. John Wiley & Sons, Inc., Hoboken, New Jersey, pp. 163–173. Knowles, L.L., Carstens, B.C., 2007. Delimiting species without monophyletic gene trees. Syst. Biol. 56, 887–895. Knowles, L.L., Kubatko, L.S., 2010. Chapter 1 Estimating species trees: an introduction to concepts and models. In: Knowles, L.L., Kubatko, L.S. (Eds.),

J.L. Malaney et al. / Molecular Phylogenetics and Evolution 114 (2017) 137–152 Estimating Species Trees: Practical and Theoretical Aspects. John Wiley & Sons, Inc., Hoboken, New Jersey. Kramer-Schadt, S., Niedballa, J., Pilgrim, J.D., Schröder, B., Lindenborn, J., Reinfelder, V., Stillfried, M., Heckmann, I., Scharf, A.K., Augeri, D.M., 2013. The importance of correcting for sampling bias in MaxEnt species distribution models. Divers. Distrib. 19, 1366–1379. Krutzsch, P.H., 1954. North American jumping mice (genus Zapus). University of Kansas Publications. Museum Nat. History 7, 349–472. Kubatko, L., Gibbs, H.L., Bloomquist, E.W., 2011. Inferring species-level phylogenies and taxonomic distinctiveness using multilocus data in Sistrurus rattlesnakes. Syst. Biol. 60, 393–409. Kubatko, L.S., Gibbs, H.L., 2010. Chapter 12 Estimating species relationships and taxon distinctiveness inSistrurus Rattlesnakes using multilocus data. In: Knowles, L.L., Kubatko, L.S. (Eds.), Estimating Species Trees Practical and Theoretical Aspects. John Wiley & Sons, Inc., Hoboken, New Jersey. Lackey, R.T., 2007. Science, scientists, and policy advocacy. Conserv. Biol. 21, 12–17. Leache, A.D., Fujita, M.K., 2010. Bayesian species delimitation in west African forest geckos (Hemidactylus fasciatus). Proc. Roy. Soc. B-Biol. Sci. 277, 3071–3077. Leache, A.D., Koo, M.S., Spencer, C.L., Papenfuss, T.J., Fisher, R.N., McGuire, J.A., 2009. Quantifying ecological, morphological, and genetic variation to delimit species in the coast horned lizard species complex (Phrynosoma). Proc. Natl. Acad. Sci. U.S.A. 106, 12418–12423. Lessa, E.P., Cook, J.A., Patton, J.L., 2003. Genetic footprints of demographic expansion in North America, but not Amazonia, during the Late Quaternary. Proc. Natl. Acad. Sci. U.S.A. 100, 10331–10334. Liu, C.R., Berry, P.M., Dawson, T.P., Pearson, R.G., 2005. Selecting thresholds of occurrence in the prediction of species distributions. Ecography 28, 385–393. Mace, G.M., 2004. The role of taxonomy in species conservation. Philos. Trans. Roy. Soc. B-Biol. Sci. 359, 711–719. Maddison, W.P., 1997. Gene trees in species trees. Syst. Biol. 46, 523–536. Maechler, M., Rousseeuw, P., Struyf, A., Hubert, M., Hornik, K., 2016. Cluster: Cluster Analysis Basics and Extensions. R package version 2.0.5. Malaney, J.L., 2012. Historical Biogeography of Boreal Mammals. Department of Biology and Museum of Southwestern Biology. University of New Mexico, Albuquerque, NM. Malaney, J.L., Conroy, C.J., Moffitt, L.A., Spoonhunter, H.D., Patton, J.L., Cook, J.A., 2013. Phylogeography of the western jumping mouse (Zapus princeps) detects deep and persistent allopatry with expansion. J. Mammal. 94, 1016–1029. Malaney, J.L., Cook, J.A., 2013. Using biogeographic history to inform conservation: the case of Preble’s meadow jumping mouse. Mol. Ecol. 22, 6000–6017. Malaney, J.L., Frey, J.K., Cook, J.A., 2012. The biogeographic legacy of an imperilled taxon provides a foundation for assessing lineage diversification, demography, and conservation genetics. Divers. Distrib. 18, 689–703. Mann, M.E., Gleick, P.H., 2015. Climate change and California drought in the 21st century. Proc. Natl. Acad. Sci. 112, 3858–3859. Mayr, E., 1956. Geographical character gradients and climatic adaptation. Evolution 10, 105–108. Mayr, E., 1963. Animal Species and Evolution. Belknap Press of Harvard University Press, Cambridge, Massachusetts. McCormack, J.E., Zellmer, A.J., Knowles, L.L., 2010. Does niche divergence accompany allopatric divergence in Aphelocoma jays as predicted under ecological speciation?: Insights from tests with niche models. Evolution 64, 1231–1244. McGarigal, K., Cushman, S., Stafford, S., 2000. Multivariate Statistics for Wildlife and Ecology Research. Springer-Verlag Inc, New York, NY, USA. Merow, C., Smith, M.J., Silander, J.A., 2013. A practical guide to MaxEnt for modeling species’ distributions: what it does, and why inputs and settings matter. Ecography 36, 1058–1069. Merriam, C.H., 1897. Three new Jumping Mice (Zapus) from the Northwest. Proc. Biol. Soc. Washington 11, 103–104. Miller, G.S., 1893. A Jumping Mouse (Zapus insignis, Miller) new to the United States. P. Biol. Soc. Washington viii, 1–8. Miller, G.S., 1911. A new Jumping-Mouse from New Mexico. Washington D.C. Proc. Biol. Soc. 24 (253–254). Morales, A.E., Jackson, N.D., Dewey, T.A., O’Meara, B.C., Carstens, B.C., 2016. Speciation with gene flow in North American Myotis bats. Syst. Biol., syw100 Moritz, C., Potter, S., 2013. The importance of an evolutionary perspective in conservation policy planning. Mol. Ecol. 22, 5969–5971. Nosil, P., 2008. Speciation with gene flow could be common. Mol. Ecol. 17, 2103– 2106. O’Meara, B.C., 2009. New heuristic methods for joint species delimitation and species tree inference. Syst. Biol., syp077 Oksanen, J., Kindt, R., Legendre, P., 2008. Vegan: Community Ecology Package. R package version 1.15-1. Olave, M., Solà, E., Knowles, L.L., 2014. Upstream analyses create problems with DNA-based species delimitation. Syst. Biol. 63, 263–271. Padial, J.M., Miralles, A., De la Riva, I., Vences, M., 2010. The integrative future of taxonomy. Front. Zool. 7, 1–14. Papadopoulou, A., Bergsten, J., Fujisawa, T., Monaghan, M.T., Barraclough, T.G., Vogler, A.P., 2008. Speciation and DNA barcodes: testing the effects of dispersal on the formation of discrete sequence clusters. Philos. Trans. Roy. Soc. B: Biol. Sci. 363, 2987–2996. Paradis, E., Claude, J., Strimmer, K., 2004. APE: analyses of phylogenetics and evolution in R language. Bioinformatics 20, 289–290. Patton, J.L., Smith, M.F., 1994. Paraphyly, polyphyly, and the nature of species boundaries in pocket gophers (Genus Thomomys). Syst. Biol. 43, 11–26.

151

Pearson, R.G., Raxworthy, C.J., Nakamura, M., Peterson, A.T., 2007. Predicting species distributions from small numbers of occurrence records: a test case using cryptic geckos in Madagascar. J. Biogeogr. 34, 102–117. Peterson, A.T., Papes, M., Soberón, J., 2008. Rethinking receiver operating characteristic analysis applications in ecological niche modeling. Ecol. Model. 213, 63–72. Phillips, S.J., Anderson, R.P., Schapire, R.E., 2006. Maximum entropy modeling of species geographic distributions. Ecol. Model. 190, 231–259. Phillips, S.J., Dudik, M., Elith, J., Graham, C.H., Lehmann, A., Leathwick, J., Ferrier, S., 2009. Sample selection bias and presence-only distribution models: implications for background and pseudo-absence data. Ecol. Appl. 19, 181–197. Pigot, A.L., Phillimore, A.B., Owens, I.P., Orme, C.D.L., 2010. The shape and temporal dynamics of phylogenetic trees arising from geographic speciation. Syst. Biol. 59, 660–673. Pons, J., Barraclough, T.G., Gomez-Zurita, J., Cardoso, A., Duran, D.P., Hazell, S., Kamoun, S., Sumlin, W.D., Vogler, A.P., 2006. Sequence-based species delimitation for the DNA taxonomy of undescribed insects. Syst. Biol. 55, 595–609. Posada, D., 2008. JModelTest: phylogenetic model averaging. Mol. Biol. Evol. 25, 1253–1256. Posada, D., Buckley, T.R., 2004. Model selection and model averaging in phylogenetics: Advantages of akaike information criterion and Bayesian approaches over likelihood ratio tests. Syst. Biol. 53, 793–808. Posada, D., Crandall, K.A., 2001. Intraspecific gene genealogies: trees grafting into networks. Trends Ecol. Evol. 16, 37–45. Preble, E.A., 1899. Revision of the jumping mice of the genus Zapus. North American Fauna 15, 1–39. Pritchard, J.K., Stephens, M., Donnelly, P., 2000. Inference of population structure using multilocus genotype data. Genetics 155, 945–959. Pritchard, J.K., Wen, X., Falusch, D., 2010. Documentation for STRUCTURE software: Version 2.3. Pyron, R.A., Burbrink, F.T., 2009. Systematics of the Common Kingsnake (Lampropeltis getula; Serpentes: Colubridae) and the burden of heritage in taxonomy. Zootaxa 2241, 22–32. Pyron, R.A., Costa, G.C., Patten, M.A., Burbrink, F.T., 2015. Phylogenetic niche conservatism and the evolutionary basis of ecological speciation. Biol. Rev. 90, 1248–1262. R Core Team, 2014. R: a language and environment for statistical computing.. Open access available at: http://cran.r-project.org or http://www.R-project.org. R Foundation for Statistical Computing, Vienna, Austria. Rambaut, A., Drummond, A., 2009. FigTree v1. 3.1. Ramey, R.R., Liu, H.P., Epps, C.W., Carpenter, L.M., Wehausen, J.D., 2005. Genetic relatedness of the Preble’s meadow jumping mouse (Zapus hudsonius preblei) to nearby subspecies of Z. hudsonius as inferred from variation in cranial morphology, mitochondrial DNA and microsatellite DNA: implications for taxonomy and conservation. Anim. Conserv. 8, 329–346. Ramey, R.R., Wehausen, J.D., Liu, H.P., Epps, C.W., Carpenter, L.M., 2006. Response to Vignieri et al. (2006): Should hypothesis testing or selective post hoc interpretation of results guide the allocation of conservation effort? Anim. Conserv. 9, 244–247. Rannala, B., 2015. The art and science of species delimitation. Curr. Zool., 61 Reddy, S., Davalos, L.M., 2003. Geographical sampling bias and its implications for conservation priorities in Africa. J. Biogeogr. 30, 1719–1727. Reid, N., Carstens, B., 2012. Phylogenetic estimation error can decrease the accuracy of species delimitation: a Bayesian implementation of the general mixed Yulecoalescent model. BMC Evol. Biol. 12, 196. Rhoads, 1894. [Title unknown.]. P. Ac. Philad., Unpaginated. Richmond, O.M.W., McEntee, J.P., Hijmans, R.J., Brashares, J.S., 2010. Is the climate right for pleistocene rewilding? Using species distribution models to extrapolate climatic suitability for mammals across continents. PLoS ONE 5, e12899. Riddle, B.R., Hafner, D.J., 1999. Species as units of analysis in ecology and biogeography: time to take the blinders off. Glob. Ecol. Biogeogr. 8, 433–441. Rissler, L.J., Apodaca, J.J., 2007. Adding more ecology into species delimitation: Ecological niche models and phylogeography help define cryptic species in the black salamander (Aneides flavipunctatus). Syst. Biol. 56, 924–942. Ronquist, F., Huelsenbeck, J., Teslenko, M., 2011. Draft MrBayes version 3.2 manual: tutorials and model summaries. Rosell, J.A., Olson, M.E., Weeks, A., De-Nova, J.A., Lemos, R.M., Camacho, J.P.R., Feria, T.P., GÃ3mez-Bermejo, R., Montero, J.C., Eguiarte, L.E., 2010. Diversification in species complexes: Tests of species origin and delimitation in the Bursera simaruba clade of tropical trees (Burseraceae). Mol. Phylogenet. Evol. 57, 798– 811. Ryder, O.A., 1986. Species conservation and systematics: the dilemma of subspecies. Trends Ecol. Evol. 1, 9–10. Satler, J.D., Carstens, B.C., Hedin, M., 2013. Multilocus species delimitation in a complex of morphologically conserved trapdoor spiders (Mygalomorphae, Antrodiaetidae, Aliatypus). Syst. Biol. 62, 805–823. Schoener, T.W., 1968. Anolis lizards of Bimini: resource partitioning in a complex fauna. Ecology 49, 704–726. Scott, J.M., Rachlow, J.L., Lackey, R.T., Pidgorna, A.B., Aycrigg, J.L., Feldman, G.R., Svancara, L.K., Rupp, D.A., Stanish, D.I., Steinhorst, R.K., 2007. Policy advocacy in science: Prevalence, perspectives, and implications for conservation biologists. Conserv. Biol. 21, 29–35. Shafer, A., Wolf, J.B., 2013. Widespread evidence for incipient ecological speciation: a meta-analysis of isolation-by-ecology. Ecol. Lett. 16, 940–950.

152

J.L. Malaney et al. / Molecular Phylogenetics and Evolution 114 (2017) 137–152

Simpson, G.G., 1951. The species concept. Evolution 5, 285–298. Sites, J.W., Marshall, J.C., 2004. Operational criteria for delimiting species. Annu. Rev. Ecol. Evol. Syst. 35, 199–227. Soberón, J., 2007. Grinnellian and Eltonian niches and geographic distributions of species. Ecol. Lett. 10, 1115–1123. Sukumaran, J., Knowles, L.L., 2017. Multispecies coalescent delimits structure, not species. Proceedings of the National Academy of Sciences, 201607921. Sullivan, J., Demboski, J., Bell, K., Hird, S., Sarver, B., Reid, N., Good, J., 2014. Divergence with gene flow within the recent chipmunk radiation (Tamias). Heredity. Tamura, K., Peterson, D., Peterson, N., Stecher, G., Nei, M., Kumar, S., 2011. MEGA5: Molecualr Evolutionary Genetics Analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol. Biol. Evol. 28, 2731–2739. Templeton, A.R., Crandall, K.A., Sing, C.F., 1992. A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping and DNA sequence data. III. Cladogram estimation. Genetics 132, 619–633. Toprak, Z., Pfeil, B.E., Jones, G., Marcussen, T., Ertekin, A.S., Oxelman, B., 2016. Species delimitation without prior knowledge: DISSECT reveals extensive cryptic speciation in the Silene aegyptiaca complex (Caryophyllaceae). Mol. Phylogenet. Evol. 102, 1–8. USFWS, 2010. Endangered and Threatened Wildlife and Plants; Revised Critical Habitat for the Preble’s Meadow Jumping Mouse in Colorado; Final Rule. In: Interior, D.o.t. (Ed.). Federal Register, pp. 78428–78483. USFWS, 2014. Endangered and Threatened Wildlife and Plants; Determination of Endangered Status for the New Mexico Meadow Jumping Mouse Throughout Its Range. In: Interior, D.o.t. (Ed.). Federal Register, pp. 33119–33137. Veloz, S.D., 2009. Spatially autocorrelated sampling falsely inflates measures of accuracy for presence-only niche models. J. Biogeogr. 36, 2290–2299. Venables, W.N., Ripley, B.D., 2002. Modern Applied Statistics With S. Springer, New York. Vignieri, S.N., Hallerman, E.M., Bergstrom, B.J., Hafner, D.J., Martin, A.P., Devers, P., Grobler, P., Hitt, N., 2006. Mistaken view of taxonomic validity undermines

conservation of an evolutionarily distinct mouse: a response to Ramey et al. (2005). Animal Conservation 9, 237–243. Warren, D.L., Glor, R.E., Turelli, M., 2008. Environmental niche equivalency versus conservatism: quantitative approaches to niche evolution. Evolution 62, 2868– 2883. Warren, D.L., Glor, R.E., Turelli, M., 2010. ENMTools: a toolbox for comparative studies of environmental niche models. Ecography 33, 607–611. Warren, D.L., Seifert, S.N., 2011. Ecological niche modeling in Maxent: the importance of model complexity and the performance of model selection criteria. Ecol. Appl. 21, 335–342. Welch, J.J., Bromham, L., 2005. Molecular dating when rates vary. Trends Ecol. Evol. 20, 320–327. Wiens, J.J., 2007. Species delimitation: New approaches for discovering diversity. Syst. Biol. 56, 875–878. Wiens, J.J., Ackerly, D.D., Allen, A.P., Anacker, B.L., Buckley, L.B., Cornell, H.V., Damschen, E.I., Jonathan Davies, T., Grytnes, J.A., Harrison, S.P., 2010. Niche conservatism as an emerging principle in ecology and conservation biology. Ecol. Lett. 13, 1310–1324. Wiens, J.J., Graham, C.H., 2005. Niche conservatism: Integrating evolution, ecology, and conservation biology. Annu. Rev. Ecol. Evol. Syst. 36, 519–539. Wiley, E.O., Lieberman, B.S., 2011. Phylogenetics: Theory and Practive of Phylogenetic Systematics. John Wiley & Sons Inc, Hoboken, New Jersey. Yang, Z., Rannala, B., 2014. Unguided species delimitation using DNA sequence data from multiple loci. Mol. Biol. Evol., msu279 Yang, Z.H., Rannala, B., 2010. Bayesian species delimitation using multilocus sequence data. Proc. Natl. Acad. Sci. U.S.A. 107, 9264–9269. Zhang, C., Zhang, D.X., Zhu, T., Yang, Z., 2011. Evaluation of a Bayesian coalescent method of species delimitation. Syst. Biol. 60, 747–761. Zhang, J., Kapli, P., Pavlidis, P., Stamatakis, A., 2013. A general species delimitation method with applications to phylogenetic placements. Bioinformatics 29, 2869–2876. Zimmermann, E.A.W., 1780. Geographische Geschichte der Menschen und der Vierfussiger Thiere. p. 432.