Specimens examined and genes sequenced

Specimens used in this study were wild-caught following guidelines approved by the American Society of Mammalogists (Sikes et al. 2011) or obtained via tissue loans and represent localities sampled across the known distribution of H. chapmani and species representing other members of the H. alfaroi species group (H. alfaroi, H. melanotis, H. rostratus, and H. saturatior; Fig. 1). Taxonomy follows Musser and Carleton (2005) with nomenclatural updates from Weksler et al. (2006). A total of 79 individuals was used in this study, of which 72 and 39 individuals were sequenced for Cytb and intron Fgb-I7, respectively. In addition, 7 Cytb sequences were obtained from GenBank (Appendix I).

DNA extraction, amplification, and sequencing

Total genomic DNA was extracted from liver tissue frozen or preserved in 95% ethanol either following Fetzner (1999) phenol–chloroform method, or using the QIAGEN DNeasy tissue kit (cat. no. 69504; Qiagen, Valencia, California). Amplification of Cytb and Fgb-I7 was performed via polymerase chain reaction (PCR) with negative controls used for all amplifications. The complete Cytb gene was amplified with the primers MVZ-05 and MVZ-14-M (modified from Smith and Patton 1993 by Arellano et al. 2005) and internal primers MVZ-45, MVZ-16 (Smith and Patton 1993). Fgb-I7 was amplified with primers B17 and Bfib (Wickliffe et al. 2003). For Cytb, the PCR master mix contained 1.0 μl of template DNA, 1.0 μl of deoxynucleotide triphosphates (dNTPs; 1.25 mM), 0.5 μl of each primer (100 μM), 3.0 μl of MgCl₂ (25 mM), 11.85 μl of distilled H₂0, and 0.15 μl of Taq polymerase. For Fgb-I7, reactions included 3.0 μl of template DNA, 1.7 μl of dNTPs (1.25 mM), 2.5 μl of each primer (100 μM), 1.7 μl of MgCl₂ (25 mM), 0.8 μl of GeneAmp 10× PCR buffer, 13.7 μl of high-performance liquid chromatograpy-H₂O, and 0.125 μl of Platinum Taq polymerase (Promega Corp., Madison, Wisconsin). Thermal profiles for Cytb were: 3 min at 94°C, 39 cycles of 1 min at 94°C, 1 min at 50°C, 1 min at 72°C, and 5 min at 72°C followed by a soak at 4°C. For Fgb-I7, a hot start of 15 min at 85°C was used before the addition of dNTPs; this was followed by 10 min at 94°C, 32 cycles of 1 min at 94°C, 1 min at 65°C, and 1 min at 72°C, and a soak at 4°C. PCR products were purified either with a Gene-Clean PCR purification kit (Bio 101, La Jolla, California) or by using a Millipore (Billerica, Massachusetts) multiscreen PCR 96-well filtration system (cat. no. MANU03050). Sequencing reactions of purified PCR products were done with the Perkin–Elmer ABI PRISM dye terminator cycle sequencing ready reaction kit (Applied Biosystems, Foster City, California). Excess dye terminator was removed using a Sephadex 50G solution (3 g/50 ml H₂O) or with a Millipore multiscreen filter plate (cat. no. MAHVN4510). Light- and heavy-strand sequences were determined with an ABI 3100 automated sequencer (Applied Biosystems) housed in the DNA Sequencing Center at Brigham Young University or by Macrogen Inc., Seoul, Korea ( http://www.macrogen.com). Sequences were edited manually using Sequencher version 4.1.1 and 4.1.2 (Gene Codes Corp., Ann Arbor, Michigan).

Phylogenetic analyses of the Cytb data set

Alignment for Cytb was done by translating nucleotide sequences into amino acids with Codon Code Aligner v2.0.6 (Codon Code Corp., Dedham, Massachusetts). Unique haplotypes were identified with TCS v1.21 (Clement et al. 2000) and models of nucleotide substitution and genetic variation parameters that best fit our data were selected using jModelTest v1.1 (Posada 2008 using the Akaike information criterion). The model of evolution selected was TVM + Γ (Posada and Crandall 1998). Base frequencies were A = 0.3332, C = 0.3209, G = 0.0981, and T = 0.2480; transversion rates were (A-C) 0.2990, (A-G) 2.4623, (A-T) 0.3909, (C-G) 0.2127, (C-T) 2.4623, and the gamma shape parameter (Γ) was 0.2310. Maximum likelihood (ML—Felsenstein 1981) and Bayesian inference (BI—Yang and Rannala 1997) optimality criteria were used to estimate relationships among taxa using RAxML v7.4.8 (Stamatakis 2006) and MrBayes 3.1.2 (Ronquist and Huelsenbeck 2003), respectively.

Handleyomys alfaroi was designated as the outgroup for our phylogenetic analyses following its current taxonomic position as the sister lineage to H. chapmani and H. rostratus (Weksler et al. 2006). We allowed H. rostratus and H. melanotis (its presumed sister lineage) to be part of the ingroup along with H. saturatior; this as an alternative test of monophyly for H. chapmani (Nixon and Carpenter 1993).

For BI, 2 analyses with 3 chains were run independently for 10 million Metropolis coupled Markov chain Monte Carlo (MCMC) generations using the default priors on model parameters starting from a random tree. For all analyses, a tree was sampled every 2,000 generations. Stationarity was determined by monitoring the fluctuating value of the likelihood parameters using Tracer v1.4 (Rambaut and Drummond 2007). All the trees before stationarity were discarded as “burn in.” For ML, a heuristic search starting from a random tree was conducted with 1,000 replicates using RAxML v7.4.8 (Stamatakis 2006). Kimura 2-parameter (Kimura 1980) genetic distances were calculated to assess within- and among-species genetic divergence using PAUP v4.0b10 (Swofford 2002) as they are directly comparable with distance values reported in treatments dealing with phylogeny reconstruction or species definitions of mammals (Smith and Patton 1993; Baker and Bradley 2006; Tobe et al. 2010).

Phylogenetic analyses of the Fgb-I7 data set

Alignment of Fgb-I7 data was done using the software MUSCLE (Edgar 2004). The model of evolution selected by jModelTest v1.1 as most appropriate for Fgb-I7 was HKY (Hasegawa et al. 1985). Base frequencies were A = 0.2949, C = 0.1725, G = 0.2035, and T = 0.3290; transition/transversion ratio = 1.2408. Phylogeny estimation was done as for Cytb data set for ML and BI.

Phylogenetic analyses of the combined data set

Before combining the data partitions, we assessed the level of disagreement between the Cytb and Fgb-I7 data sets with the incongruence length difference test (ILD—Farris et al. 1995; see also Hipp et al. 2004) using simple taxon addition, nearest-neighbor interchange branch swapping, and a heuristic search using 1,000 replicates in PAUP v4.0b10 (Swofford 2002). Initially, the test was run by comparing 1,143 base pairs (bp) of Cytb with 621 bp of Fgb-I7, which resulted in rejecting the null hypothesis of data homogeneity (P = 0.01). Then, Cytb was reduced to its first 621 bp and to its last 621 bp to match the length of Fgb-I7. In both cases, these tests failed to reject the null hypothesis of data homogeneity (P = 0.09). This inconsistency highlights some of the criticisms of the ILD test (Yoder et al. 2001; Barker and Lutzoni 2002; Hipp et al. 2004). Alternatively, studies have demonstrated that total evidence may provide better resolution than separate analyses that are not fully resolved (Chippindale and Wiens 1994; Jackman et al. 1997), especially when the conflict is small and most regions of the tree are shared between partitions (Wiens 1998). Therefore, we followed Wiens (1998) methodology for data combinability and analyzed each partition separately to identify parts of the tree where there was incongruence; then combined the data sets and considered the conflicted branches with caution. All major haploclades recovered by Cytb were represented in the combined data set. For BI and ML combined analysis, the partition substitution models formerly selected were specified. Combined data analyses were run with the same settings as described for the Cytb.

Nodal support

ML branch support was determined with 1,000 nonparametric bootstrap pseudoreplicates (Felsenstein 1985). Bootstrapping was synchronized with phylogeny reconstruction in RAxML v7.4.8 (Stamatakis 2006). Clades with bootstrap proportions (BP) above 70% were considered relatively well supported (Hillis and Bull 1993). For Bayesian analyses, the posterior probabilities (pP) for individual clades were obtained by constructing a majority-rule consensus of the trees not discarded as burn-in using MrBayes 3.1.2 (Ronquist and Huelsenbeck 2003).

Topology tests

Statistical support for the tree topologies was estimated as the pP for the subset of possible trees in agreement with the topology we recovered (Huelsenbeck and Rannala 2004). Specific ML topology tests were performed with the approximately unbiased (AU) test (Shimodaira 2002) with 10 sets of 10,000 bootstrap replicates. Both tests were performed on the combined data set and run in Consel (Shimodaira and Hasegawa 2001) using site log-likelihoods calculated with PAUP v4.0b10 (Swofford 2002). The model of nucleotide substitution for the AU test was GTR with optimized parameters using RAxML.

Divergence times estimates

Phylogeny dating was assessed with the coalescent Bayesian approach for multilocus data implemented in BEAST v1.7.4 (Drummond and Rambaut 2007). For each partition, parameters for the model of nucleotide substitution were the same used for the phylogenetic analyses (unlinked substitution models). Two MCMC analyses were run for 10,000,000 generations with trees sampled every 1,000 generations. Stationary, appropriate effective sample size (ESS) and convergence of independent MCMC was visualized with Tracer v1.5. The first 1,000 trees of each run (10% respectively) were discarded as burn-in and the remaining trees were then combined to build a maximum credibility tree in TreeAnnotator v1.6. Analyses were done under the assumption of a relaxed molecular clock to account for heterogeneity of substitution rates among lineages (Arbogast et al. 2002). Using a Yule tree prior on the net rate of speciation, rates among lineages were assumed to be uncorrelated, and the rate for each branch was independently drawn from a lognormal distribution (uncorrelated lognormal model—Drummond et al. 2006). As calibration points we used fossil records for H. alfaroi and H. rostratus (= H. melanotis) from the Rancholabrean 0.3 million years ago (mya—Ferrusquía-Villafranca et al. 2010) and H. alfaroi from the late Quaternary (0.5–1.0 mya—Arroyo-Cabrales et al. 2002). This information was incorporated in the H. rostratus and H. alfaroi nodes to set hard lower bounds of 0.3 mya and 0.5 mya, respectively for the time of most recent common ancestor (tMRCA).

Delimiting species boundaries

Although species delimitation is an inherent practice in phylogenetics, until recently, implementation of species boundaries had lacked a theoretical framework on which such limits could be tested explicitly (Sites and Marshall 2003, 2004; Wiens 2007; Rogers and Gonzalez 2010). This is attributable, at least in part, to the natural subjectivity of species concepts and incompatibilities among them (de Queiroz 2007). Nevertheless, this topic is receiving more attention and hypothesis-testing methods for delimitation of species boundaries have been developed (see Wiens 2007; Camargo et al. 2012; Fujita et al. 2012).

The amount and direction of gene flow was estimated using MIGRATE-N v.3.3 (Beerli and Felsenstein 2001) as the mutation scaled effective migration rate (M) to account for the autosomal inheritance of Fgb-I7. M in turn was multiplied by the estimated effective population size (theta = Θ) to obtain the effective number of migrants per generation (Nm). F_ST estimates were used as starting values to run 3 replicate chains with 100,000 genealogies. If migration between lineages was not perceived, the degree of exclusive ancestry was quantified with the genealogical sorting index (GSI—Cummings et al. 2008). GSI ranges from 0 to 1, where values < 1 basically reflect additional coalescent events from the minimum required to unite all members of the group through a most recent common ancestor. Statistical significance for the GSI values (probability of finding that degree of exclusive ancestry in our groupings by chance) is estimated with a permutation test (Cummings et al. 2008). For BI trees (Cytb, Fgb-I7, and concatenated), GSI values were calculated for the last 100 trees of the MCMC search. For ML trees (Cytb, Fgb-I7, and concatenated), GSI values were calculated on the best tree found during the heuristic search (with RAxML v7.4.8; see “Materials and Methods”). We also calculated the GSI for the Cytb and Fgb-I7 gene topologies ensemble (GSI_T).

Geographic association of the recovered lineages was assessed with the nested clade phylogeographical analysis (NCPA—Templeton 1998) and GeoDis (Posada et al. 2000) run in their automated form as implemented in ANeCA (Panchal 2007). Although the NCPA has been criticized, most of these arguments are based on the lack of statistical assessment of uncertainty and related to inconsistencies of the inferences of complex phylogeographic histories, particularly involving high migration rates (Beaumont and Panchal 2008; see Beaumont et al. 2010 for a detailed review). However, the performance of this method has been defended (Templeton 2008, 2010a, 2010b). For the purpose of this paper, migration rates were explicitly estimated previous to this test, and under those circumstances the test provides a concrete way to describe the distribution of genetic variation over geography.

Finally, using the sequence data from both markers, we estimated the posterior probability (BpP) for a model of speciation using those clades that were characterized by a lack of migration and suggested as significantly exclusive on the basis of results of the GSI tests. These analyses were implemented in the software Bayesian phylogenetics and phylogeography (BPP v.2.0—Yang and Rannala 2010). The coalescent species delimitation method used in BPP relies on a reversible-jump Markov chain Monte Carlo (rjMCMC) for taking into account uncertainty due to unknown gene trees and ancestral coalescent processes. An equal prior probability was assigned to all species delimitation models (1–4 species, 5 species, and 6 species), and to ensure convergence of the estimates, rjMCMC was run with algorithm 0 and 1 (with fine-tune parameter ϵ = 15) starting from 3 different trees (fully resolved; 6 species, 5 species, and 1 species). The rjMCMC was run for 500,000 generations with a sampling frequency of 5 after a burn-in period of 10,000. The mean value for the ancestral population size (Θ) was estimated with MIGRATE-N v.3.3 to set a gamma prior (α, β) of Θ = (5.0, 100) and root age (tau = τ) τ0 = (2, 1,000). The Kimura 2-parameter (Kimura 1980) genetic distances for Cytb were then used for sequence divergence comparisons within and among the lineages identified.

Phylogenetic analysis of individual genes

Of the 1,143 Cytb nucleotides, 329 were variable. Both ML and BI phylogenetic analyses converged on basically the same tree topology (Fig. 2). Nodal support was high for all species-level taxa and geographically exclusive clades. A total of 65 haplotypes in the Cytb data set were identified, of which 49 represented samples of H. chapmani and 16 belonged to other species of Handleyomys. With only 3 exceptions, H. chapmani haplotypes also were exclusive by locality. These exceptions included 1 haplotype found at localities 11 (Colección de Mamíferos del Centro de Investigación en Biodiversidad y Conservación, Universidad Autónoma del Estado de Morelos [CMC] 779, CMC 782) and 7 (CMC 1052; Fig. 1), and 1 haplotype present in localities 9 (CMC 1450), 6 (Monte L. Bean Life Science Museum, Brigham Young University [BYU] 15803), and 7 (CMC 1049, CMC 1054); all within SMO. Similarly, within the OH a haplotype from locality 16 (CMC 114) was found at locality 14 (CMC 1347). Haplotypes representing H. alfaroi, H. melanotis, and H. rostratus also were exclusive by locality, except for a H. melanotis haplotype that was present at locality 33 (ASNHC 3418) and locality 34 (ASNHC 3419). Topologies generated under both ML and BI optimality criteria showed that H. chapmani is not a monophyletic group. Samples of this species were recovered in 2 divergent clades (SMO-OH and SMS; Fig. 2) with high nodal support (pP = 1.0, BP = 100 for both). Furthermore, H. saturatior was placed as the sister group to the H. chapmani SMO-OH clade (pP = 0.92, BP = 91). H. melanotis and H. rostratus were recovered as sister taxa (pP = 1.00, BP = 100), and samples of each species constituted strongly supported monophyletic assemblages (pP = 1.00, BP = 100). H. rostratus was recovered in 2 well-supported clades (pP = 1.0, BP = 91–94) corresponding to samples from east and west of the Isthmus of Tehuantepec (Fig. 2).

The Fgb-I7 data set consisted of 621 characters, of which 52 were variable. Three indels were inferred for our Fgb-I7 sequences based on H. alfaroi as the outgroup. One was assigned at position 283 (single bp deletion for H. chapmani and H. saturatior), a 2nd gap was identified at positions 382–383 (an insertion inferred for all ingroup taxa), and a 3rd indel was set at positions 407–417 (a deletion inferred for H. rostratus). There were 14 Fgb-I7 haplotypes identified by TCS (Phylogenetic network estimation using statistical parsimony; Clement et al. 2000), 8 of which were present only in H. chapmani. Of these 8, 3 were unique haplotypes and 5 were shared but exclusive by regions (SMS, SMO-OH). H. saturatior was represented by 2 haplotypes corresponding to different localities, H. alfaroi by the same haplotype found at 2 localities, H. melanotis by 1 haplotype present in all 5 localities, and H. rostratus by 2 haplotypes from a single locality.

Phylogenetic analyses of Fgb-I7 based on ML and BI optimality criteria estimated genealogies that were highly concordant, albeit less resolved than those recovered with Cytb (Fig. 3a). Samples of H. chapmani and H. saturatior grouped together as a well-supported monophyletic assemblage (pP = 1.00, BP = 100), although this clade resulted in an unresolved polytomy. Nonetheless, within this clade samples were arranged following a geographic pattern by mountain range. Samples of H. chapmani from the SMS formed 2 clades, 1 comprising CMC 1655 and CMC 1657 from El Tejocote, Guerrero (locality 20; pP = 1.00, BP = 89) and the other comprising the remaining samples. Similarly, all H. saturatior samples except ECOSCM 1231 (locality 24) also were recovered in a well-supported clade (pP = 1.00, BP = 98). H. melanotis and H. rostratus were recovered as monophyletic clades (pP = 1.00, BP = 100, for both clades). However, H. melanotis was placed as sister group to the H. chapmani and H. saturatior clade (pP = 0.90, BP = 90).

Phylogenetic analysis of combined data set

Trees estimated from the combined data set (Cytb and Fgb-I7) converged on basically the same tree topology for both ML and BI optimality criteria, as recovered by the Cytb tree (Fig. 2). Fig. 3b depicts the ML tree (lnL = −5754.025). H. chapmani was recovered as 2 polyphyletic clades (SMO-OH and SMS). Each of the clades recovered was strongly supported by Bayesian pP and ML bootstrap values (pP = 1.00, BP = 100). Additionally, the H. chapmani SMO-OH clade was placed as the sister group to H. saturatior (pP = 0.90, BP = 83).

Topology tests

The pP value for a topology that constrained clades representing H. chapmani to be monophyletic was pP = 0.142, whereas the probability of a topology with H. chapmani paraphyletic (as recovered in this study) was pP = 0.858. With the AU test, the hypothesis of H. chapmani monophyly was rejected (AU = 0.0451; P < 0.05). Log-likelihood of the constrained topology (H. chapmani monophyletic) was lnL = −5,762.05701, whereas the unconstrained topology (H. chapmani paraphyletic) was lnL = −5,754.025.

Divergence times estimates

The MCMC combined runs reached ESS above 450 for all parameters. The standard deviation of the uncorrelated lognormal relaxed clock for Cytb had a mean of 0.883 and 1.982 for Fgb-I7, indicating that both were not behaving in a clock-like manner. The mean substitution rate (per site per million years) was 0.027 for Cytb and 0.009 for Fgb-I7. On a timescale, H. alfaroi was not used to root the tree because the root is implicit as the most recent common ancestor (MRCA). The root was placed with a mean age of 2.51 mya (highest posterior density interval [95% HPD] = 1.30, 3.76). A mean divergence time of 1.45 mya was estimated for the H. chapmani SMS clade (95% HPD = 0.65, 2.40), whereas the H. chapmani SMO-OH and H. saturatior split was estimated at 1.08 mya (95% HPD = 0.54, 1.86). The divergence time estimate for H. melanotis and H. rostratus was placed at 1.53 mya (95% HPD = 0.68, 2.50). When it was not constrained as the root, H. alfaroi was positioned as sister to the H. melanotis and H. rostratus clade; the tMRCA for this clade was estimated at 2.07 mya (95% HPD = 1.08, 2.97).

Inferred species boundaries

There was no evidence of gene flow between H. chapmani (SMO-OH) and H. saturatior (Nm = 0.09 [95% HDP = 0.00–0.18]), between H. chapmani (SMO-OH) and H. chapmani (SMS; Nm = 0.07 [95% HDP = 0.00 – 0.14]), or between H. chapmani (SMS) and H. saturatior (Nm = 0.09 [95% HDP = 0.00 – 0.18]). The opposite migration estimates were equivalent and are not shown.

The GSI values for H. chapmani as currently defined (SMO-OH and SMS clades labeled as H. chapmani) averaged 0.794 (min. = 0.552, max. = 0.894; Table 1). When H. chapmani was labeled according to the 2 clades recovered in this study, the GSI values were higher for each lineage (SMO-OH averaged 0.889 [min. = 0.655, max. = 1] and SMS averaged 0.862 [min. = 0.604, max. = 1]). The GSI values for H. saturatior were smaller (min. = 0.379, max. = 1, average = 0.779) than those calculated for each H. chapmani clade. For each basal lineage recovered in our study, Cytb and concatenated data topologies consistently recovered values of 1 (achieved monophyly); and weighted GSI analyses (Cytb and Fgb-I7 topologies combined; GSI_T) ranged from 0.664 to 0.827. All GSI statistics had significant P-values (< 0.0004). Fgb-I7 trees yielded lower GSI values for all the groupings (Table 1).

The species delimitation analyses (BPP) strongly supported a model of 5 speciation events (6 species; Fig. 4), corresponding to H. alfaroi, H. melanotis, H. rostratus, H. chapmani SMS, H. saturatior, and H. chapmani SMO-OH (BpP = 0.99620). Under this model, H. chapmani SMO-OH, H. chapmani SMS, and H. saturatior have a BpP = 1.00000; H. alfaroi a BpP = 0.99999, and H. rostratus and H. melanotis a BpP = 0.99666. In contrast, a model of 5 species had a BpP = 0.00374, and a model of < 5 species had a BpP = 0.00005. To examine the effect of excessive a priori subdivision, we further split H. chapmani SMO-OH to create a model with 7 species (all of the above species plus H. chapmani OH populations as the 7th). This model had a much lower probability BpP = 0.09518 than our 6-species speciation model (BpP = 0.99620).

The Cytb haplotype networks we identified corresponded to the clades recovered by ML and BI analyses. Samples of H. chapmani were represented in 2 separate networks (SMS and SMO-OH). The NCPA indicated that haplotypes of H. chapmani corresponding to the SMS and SMO-OH clades are geographically isolated genetic clusters (P = 0.0277 and P = 0.0411 respectively). Within SMS (Geodis Dc), the southern Oaxaca and western Guerrero haplotypes are the most geographically restricted (P = 0.0121). Within SMO, haplotypes found in San Luis Potosí connected to those in Hidalgo but with a significantly large nested clade distance (Dn; P = 0.0249). In contrast, haplotypes from northern Oaxaca (OH) were recovered as a separate genetic unit from the rest of H. chapmani SMO (P = 0.0010; Figs. 1 and 2).