Handleyomys chapmani (Chapman's Handley's mouse) is a Mexican endemic rodent inhabiting humid montane forest of the Sierra Madre Oriental (SMO), the Oaxacan Highlands (OH), and the Sierra Madre del Sur (SMS). The systematic status of populations currently classified as H. chapmani has been problematic and to date evolutionary relationships among populations remain unresolved. In this study we use sequences from the mitochondrial cytochrome-b gene (1,143 base pairs [bp]) and intron 7 of the beta fibrinogen gene (621 bp) to reconstruct a phylogeny, estimate divergence times, and assess patterns of sequence variation over geography among samples of H. chapmani. This species was recovered as 2 monophyletic clades corresponding to the SMO-OH and SMS mountain ranges. Moreover, H. saturatior, the purported sister taxon to H. chapmani, was consistently recovered as the sister lineage to the SMO-OH clade, rendering H. chapmani paraphyletic. The geographic distribution of the 2 H. chapmani clades and of H. saturatior strongly correlate with the geographic extent of the SMO-OH, SMS, and the Trans-Isthmian Highlands (highlands east of the Isthmus of Tehuantepec through Central America) mountain ranges. Divergence times associate their isolation to late Pleistocene climatic changes that likely were reinforced by barriers such as the Isthmus of Tehuantepec, the Tehuacán–Cuicatlán Valley, and the Central Valleys of Oaxaca. The fact that populations of H. chapmani represent 2 independent evolutionary lineages results in a substantial reduction in the distributional range for both entities. Therefore, the conservation status of H. chapmani should be re-evaluated.
Handleyomys chapmani (ratón de Handley de Chapman) es un roedor endémico de México con distribución en la Sierra Madre Oriental (SMO), Sierra Norte de Oaxaca (OH) y Sierra Madre del Sur (SMS). El estatus taxonómico de las poblaciones actualmente clasificadas como H. chapmani ha sido problemático y hasta la fecha, las relaciones evolutivas entre dichas poblaciones continúan sin resolverse. En este estudio, usamos secuencias del gen mitocondrial citocromo b (1143pb) y del intron 7 del gen beta fibrina (621pb) para estimar una filogenia del grupo, tiempos de divergencia y analizar los patrones de variación genética entre poblaciones de H. chapmani en un sentido geográfico. H. chapmani fue recuperado en 2 clados monofiléticos correspondientes a los sistemas montañosos de la SMO-OH y SMS. Además, H. saturatior (ratón de Handley de bosque nublado), reconocido como el grupo hermano de H. chapmani, fue consistentemente recuperado como el linaje hermano al clado de las SMO-OH; revelando a H. chapmani como un taxón parafilético. La distribución geográfica de los 2 clados en H. chapmani y H. saturatior muestra una fuerte correlación con la extensión geográfica de la SMO-OH, la SMS y las Tierras Altas Trans-Istmicas (TIH; tierras altas al este del Istmo de Tehuantepec en Chiapas y América Central). Los tiempos de divergencia asocian el aislamiento de éstas entidades con cambios climáticos del Pleistoceno superior, que posiblemente fue reforzado por barreras geográficas como el Istmo de Tehuantepec, el Valle Tehuacán-Cuicatlán y los Valles Centrales de Oaxaca. El hecho de que las poblaciones de H. chapmani constituyan 2 entidades evolutivas, tiene como consecuencia la reducción significativa del rango de distribución de estos 2 linajes. Por lo tanto, el estatus de conservación de H. chapmani debe ser reevaluado.
Recent studies using molecular data have demonstrated that rodent populations from different mountain ranges in Mexico exhibit considerable levels of genetic differentiation (Sullivan et al. 1997; Harris et al. 2000; Arellano et al. 2005; León-Paniagua et al. 2007; Rogers et al. 2007; Rogers and González 2010; Vallejo and González-Cózatl 2012; Hardy et al. 2013). Within the Handleyomys alfaroi group, some forms are confined to medium- and high-elevation forests in different mountain ranges of Mesoamerica (Musser and Carleton 2005). Taxonomically, this species group had been included in Oryzomys, but a systematic evaluation conducted by Weksler et al. (2006) proposed that the genus be restricted to the “palustris group” and the remaining 10 clades were elevated to generic rank, including the “alfaroi group,” which was provisionally assigned to Handleyomys. As a result, we refer to members of the H. alfaroi group as Handleyomys rather than Oryzomys throughout this paper.
The H. alfaroi group (Goldman 1918; Hall 1981), as currently defined (Weksler et al. 2006), is a complex of 6 species that includes H. alfaroi (Alfaro's Handley's mouse), H. chapmani (Chapman's Handley's mouse), H. melanotis (Black-eared Handley's mouse), H. rostratus (Long-nosed Handley's mouse), H. rhabdops (Highland Handley's mouse), and H. saturatior (Cloud Forest Handley's mouse). Although previous workers (Goldman 1915, 1918; Hall 1981; Musser and Carleton 1993, 2005) had retained H. melanotis and H. rostratus in the melanotis group, Weksler et al. (2006) included these 2 species in the H. alfaroi group. The systematics of this species group has been controversial and, over time, this complex has included from 5 (Goldman 1918) to 12 species (Allen 1891, 1913; Allen and Chapman 1897; Merriam 1901; Goldman 1915). Within this complex, the taxonomy of the Mexican endemic H. chapmani also has been problematic. The 1st specimens referable to this taxon were collected near Jalapa, Veracruz, in the Sierra Madre Oriental (SMO), and were regarded as H. melanotis by Allen and Chapman (1897). Later, Thomas (1898) referred to these specimens as H. chapmani. In 1901, Merriam recognized H. chapmani (sensu Thomas 1898) and described specimens from northern Oaxaca (Oaxacan Highlands; OH) as H. c. caudatus and those from Puebla as H. c. dilutior (SMO). Goldman (1915) then described specimens from Guerrero and southern Oaxaca in the Sierra Madre Sur (SMS) as H. guerrerensis. In his revision of North American rice rats, Goldman (1918) retained H. guerrerensis as a full species, but relegated H. chapmani and the 2 subspecies contained therein (caudatus and dilutior) as subspecies of H. alfaroi. Specimens collected by Dalquest (1951) from Tamaulipas and San Luis Potosí (SMO) were described as a new subspecies of H. alfaroi (H. a. huastecae). Interestingly, Goodwin (1969) recognized 2 different forms of Handleyomys in the mountains of northeastern Oaxaca (OH); the larger specimens were described as H. caudatus, whereas the smaller form was viewed as H. a. chapmani. Additionally, specimens of H. guerrerensis from southern Oaxaca were relegated to a subspecies of H. alfaroi (H. a. guerrerensis—Goodwin 1969). Engstrom (1984) reported a unique karyotype for H. caudatus and recognized it as distinct from H. alfaroi. More recently, Musser and Carleton (1993, 2005) considered that all described forms restricted to cloud forests in the SMO, OH, and SMS (H. a. chapmani, H. a. dilutior, H. a. guerrerensis, H. a. huastecae, and H. caudatus) were conspecific and classified under the name H. chapmani with H. saturatior as the sister group.
Given that H. chapmani is distributed allopatrically across a series of mountain ranges in northern Mesoamerica and has a complicated taxonomic history, we used deoxyribonucleic acid (DNA) sequence data from the mitochondrial gene cytochrome-b (Cytb) and the nuclear intron 7 of the beta fibrinogen (Fgb-I7) as a first approach to estimate phylogenetic relationships among populations of H. chapmani. Other members of the H. alfaroi group for which tissue samples were available (H. alfaroi, H. melanotis, H. rostratus, and H. saturatior) also were included to estimate the phylogenetic affinities relative to H. chapmani. Specifically, we use our sequence data to test the proposal of Musser and Carleton (1993, 2005) that all forms of Handleyomys restricted to cloud forest elevations of the SMO, OH, and SMS and currently considered as conspecific forms of H. chapmani (H. a. chapmani, H. a. dilutior, H. a. guerrerensis, H. a. huastecae, and H. caudatus) represent a monophyletic assemblage. Also, we test the hypothesis that H. chapmani and H. saturatior represent sister species (Musser and Carleton 1993, 2005).
Materials and Methods
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 MgCl2 (25 mM), 11.85 μl of distilled H20, 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 MgCl2 (25 mM), 0.8 μl of GeneAmp 10× PCR buffer, 13.7 μl of high-performance liquid chromatograpy-H2O, 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 H2O) 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.
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).
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). FST 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 (GSIT).
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).
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; GSIT) 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).
Our molecular phylogeny demonstrates that H. chapmani is comprised of 2 nonsister lineages that are restricted to different mountain systems (SMO-OH and SMS). Moreover, the SMO-OH clade is the sister group to H. saturatior, which occurs to the east of the Isthmus of Tehuantepec in the highlands of Chiapas and Central America (TIH). Together, these 3 lineages constitute a well-supported monophyletic assemblage. By extension, our phylogenetic analyses do not support the proposal of Musser and Carleton (1993, 2005) that all forms of H. chapmani (H. a. chapmani, H. a. dilutior, H. a. guerrerensis, H. a. huastecae, and H. caudatus) restricted to cloud forest habitat in the SMO, OH, and SMS are conspecific.
The mean Cytb genetic distance among clades of H. chapmani from the SMO-OH and SMS was 6.5%, whereas values between H. chapmani SMO-OH and H. saturatior and between H. chapmani SMS and H. saturatior were 6.0% and 6.9%, respectively. These values are comparable with those among many cryptic species of mammals (Baker and Bradley 2006). Also, the degree of genetic differentiation is in agreement with our divergence time estimates, which suggests that the 2 lineages of H. chapmani (SMO-OH and SMS) and H. saturatior are the most recently derived lineages within the H. alfaroi group (1.08–1.45 mya). According to Arroyo-Cabrales et al. (2002), Ceballos et al. (2010), and Ferrusquía-Villafranca et al. (2010), H. rostratus and H. alfaroi were well-differentiated forms in the Pleistocene fauna of Mexico. This proposal is consistent with the estimated MRCA for H. alfaroi–H. rostratus and H. melanotis (2.07 mya) and the relatively large percent Cytb sequence divergence between them (13.2%).
Although Fgb-I7 has been useful in resolving intrageneric relationships in other rodent groups (Matocq et al. 2007; Hanson and Bradley 2009), it typically has a slower substitution rate than Cytb in mammals (Wickliffe et al. 2003). We interpret the lack of resolution in the Fgb-I7 topology as a case of incomplete lineage sorting. This is supported by the lack of evidence for gene flow and the Fgb-I7 GSI values showing a substantial amount of exclusivity for each H. chapmani lineage (0.65 to 0.83) despite the partially resolved phylogeny.
Lack of detectable gene flow between H. chapmani clades SMO-OH and SMS (Nm = 0.11) and between any of these 2 lineages and H. saturatior (average Nm = 0.13) also support the notion that these groups represent separate biological species. Similarly, the geographic distributions of these 3 lineages are allopatric as supported by the NCPA analysis. The GSI values for H. chapmani clades SMO-OH and SMS showed considerable amounts of exclusive ancestry (mean GSI values of 0.889 and 0.862, respectively) and achieved monophyly (GSI = 1) with Cytb and concatenated data topologies. Moreover, the GSI values for the 2 H. chapmani clades were consistently larger than for H. saturatior, whose average GSI value was 0.779. Accordingly, BPP assigned the highest probability to a speciation model in which H. chapmani SMO-OH and H. chapmani SMS constitute 2 separate species (BpP = 1.0 for each lineage). This interpretation is also consistent with the phylogenetic species concept (Cracraft 1989). Therefore, we regard the SMS evolutionary lineage of H. chapmani as an unrecognized species on the basis of our molecular genealogy and morphological differences (Goldman 1915, 1918) from the SMO-OH clade.
Goldman (1915) described individuals from Omiltemi, Guerrero (SMS; locality 21; CMC 455) as H. guerrerensis. Later, he incorporated individuals from southern Oaxaca (SMS; ∼10 km E locality 17; CMC 943) and extended the geographic distribution of this taxon to the “forested Pacific slopes of the Sierra Madre del Sur in Guerrero and Oaxaca” (Goldman 1918:69). In comparison with H. chapmani from the SMO-OH, Goldman (1918:70) described skulls representing the SMS form as “smaller and flatter; zygomata tending to curve evenly outward, the sides less nearly parallel; sides of rostrum more tapering anteriorly; ascending branches of premaxillae usually broader posteriorly; maxillary arms of zygoma more slender; incisors smaller.” Goodwin (1956) acknowledged the morphological uniqueness of the SMS form described by Goldman (1918) and retained it as species. In a review of the mammals of Oaxaca, Goodwin (1969) included individuals from the SMS localities 17 (CMC 943), 18 (CMC 925, CMC 930, CMC 931, CMC 932, CMC 927), and ∼80 km NE (by road) locality 20 (CMC 1655, CMC 1656, CMC 1657) and compared them with H. chapmani from the SMO-OH including locality15 (CMC 113, CMC 115) and locality 16 (SMO-OH; CMC 114, CMC 117, CMC 119). Although Goodwin (1969) acknowledged the morphological features underlined by Goldman (1918), he relegated H. guerrerensis to a subspecies of H. alfaroi.
Because the name chapmani first was assigned to voucher specimens of Handleyomys from Xalapa, Veracruz (Thomas 1898—but originally described as H. melanotis by Allen and Chapman ), and our sampling included 1 specimen from this locality (CMC1450; locality 9) plus 4 more from a nearby location (CMC 1495, CMC 1497; locality 8; CMC 1353, CMC 1490; locality 10), we propose that the SMO-OH clade should retain the name chapmani. Determining a valid name for the SMS clade would require sequence data from specimens collected at or near type localities of names currently in synonymy. Goldman (1915) first used the term guerrerensis and Goodwin (1956, 1969) preserved the name guerrerensis. Therefore, taking into account that our sampling of the SMS included specimens from the type locality of guerrerensis (CMC 455; locality 21), we propose that guerrerensis is the name with priority and should be applied to the H. chapmani SMS clade. H. saturatior originally was described as a subspecies of H. chapmani (Merriam 1901) and later was retained as a subspecies of H. alfaroi (Goldman 1918). Our results support recognition by Musser and Carleton (1993, 2005) of H. saturatior as a species-level taxon.
There is general agreement that the highlands of the SMO, OH, SMS, and the different mountain ranges in the TIH represent different biogeographic provinces (Halffter 1987; Liebherr 1994; Marshall and Liebherr 2000; Contreras-Medina et al. 2007; Morrone 2010). Overall, the distributional patterns observed in these studies are supported by a variety of taxa, including plants and various animal groups (Luna-Vega et al. 2001; García-Moreno et al. 2004, 2006; Contreras-Medina et al. 2007; Puebla-Olivares et al. 2008; Bryson et al. 2011). The SMS and SMO provinces are thought to have been separated by intense volcanism in the Miocene (∼15 mya) during the formation and migration of the Mexican Transvolcanic Belt, with continuing volcanism until ∼3.5 mya (Ferrari et al. 1999). Similarly, the highlands south of the Isthmus of Tehuantepec were repeatedly isolated from mountain ranges to the north and east, with the most recent marine incursion thought to have occurred in the late Pliocene ∼3.6 mya (Maldonado-Koerdell 1964; Beard et al. 1982; Coates and Obando 1996). Climatic changes during the Pleistocene (∼2.5 mya) could have reinforced the isolating effects of a low-lying isthmus (Toledo 1982), as supported by our divergence time estimates.
Overall, levels of genetic differentiation within the H. chapmani–H. saturatior complex are in agreement with 3 main clades that occur in isolated mountain ranges (SMO-OH, SMS, and TIH). Therefore, it is reasonable to assume that these lineages have been subjected to similar historical genetic isolation and diversification, as have other montane rodent taxa such as Peromyscus (Sullivan et al. 1997; Harris et al. 2000), Reithrodontomys (Arellano et al. 2005; Hardy et al. 2013), Habromys (León-Paniagua et al. 2007; Rogers et al. 2007), and Glaucomys (Kerhoulas and Arbogast 2010), as well as a variety of other vertebrate taxa (see Almendra and Rogers 2012 for a recent summary). However, the degree of divergence of the splits among the 3 main lineages of the H. chapmani–H. saturatior complex is not completely consistent with the general patterns observed in other taxonomic groups inhabiting montane systems. Although the lowlands of the Tehuacán–Cuicatlán Valley and the Central Valleys of Oaxaca separate the SMO-OH and SMS, biogeographically, it would be more plausible to expect a closer relationship between the 2 lineages of H. chapmani (SMO-OH and SMS). This is because their distributional ranges are closer to each other than either is to H. saturatior (TIH).
It has been suggested that the Isthmus of Tehuantepec represents the deepest biogeographic break for closely related taxa of rodents with a geographic distribution along the highlands of México and Central America (Sullivan et al. 2000). Likewise, genetic differentiation recovered herein has been replicated for other rodent clades whose distributions span the Isthmus of Tehuantepec: Peromyscus (Sullivan et al. 1997); Reithrodontomys (Arellano et al. 2005; Hardy et al. 2013), Habromys (León-Paniagua et al. 2007), and Neotoma (Edwards and Bradley 2002) as well as other highland taxa (birds—Weir et al. 2008; Barber and Klicka 2010; and reptiles—Castoe et al. 2009). Nevertheless, our data show that within the H. chapmani–H. saturatior complex, the deepest split corresponds to the Tehuacán–Cuicatlán Valley and the Central Valleys of Oaxaca, rather than the Isthmus of Tehuantepec. The isthmus has played an important role in the evolutionary diversification of H. chapmani (SMO-OH)–H. saturatior, as noted by Musser and Carleton (1993, 2005).
It is interesting to note that even though the distribution area of the H. chapmani SMO-OH clade includes 2 mountain ranges that are split by the Rio Santo Domingo valley in Oaxaca, samples from each mountain system were not separated in our phylogenetic analyses. This pattern is consistent with that of other rodent species that are continuously distributed along highlands of the SMO and OH, but with no apparent genetic differentiation between samples occurring on each mountain system (i.e., the Mexican harvest mouse, Reithrodontomys mexicanus—Arellano et al. 2005). There are, however, examples of other groups of rodents in which the Rio Santo Domingo has played an important role in the diversification of populations on either side of this geological barrier (Jico deermouse, Habromys simulatus [SMO]; Chinanteco deermouse, H. chinanteco [OH]—Carleton et al. 2002; Rogers et al. 2007; Nelson's big-toothed deermouse, Megadontomys nelsoni [SMO]; Oaxacan big-toothed deermouse, M. cryophilus [OH]—Vallejo and González-Cózatl 2012). The only evidence of differentiation between samples of Handleyomys chapmani from SMO and OH was generated by our NCPA analysis, which found that haplotypes from northern Oaxaca (OH) constitute an allopatric genetic unit from the rest of H. chapmani SMO (P = 0.0010).
Mexico is considered a biodiversity hot spot for mammals, both in terms of species richness and endemism (Ceballos et al. 1998; Ceballos and Ehrlich 2006; Ceballos 2007; Giam et al. 2011). Tropical montane cloud forest is the most diverse vegetation type in Mexico, but comprises only 1% of the land surface of the country (Pedraza and Williams-Linera 2003). Unfortunately, cloud forest habitat has had a loss of 41% of its original land area, and of what remains, more than 52% is degraded (Mas et al. 2009; Sánchez Colón et al. 2009). As a result of habitat loss, H. saturatior is currently listed as near threatened (Reid et al. 2008) and H. rhabdops is listed as vulnerable (Reid and Vázquez 2008). Despite the high rates of cloud forest deforestation, H. chapmani is listed as least concern mainly because of its relatively large distribution (Castro-Arellano and Vázquez 2008). However, the fact that we recovered 2 evolutionary units within H. chapmani results in a substantial reduction in range for both lineages. As a result, the conservation status of H. chapmani should be re-evaluated.
We acknowledge financial support from the Department of Biology, the Office of Research and Creative Activities and the Monte L. Bean Life Science Museum at Brigham Young University (to DSR), and the Consejo Nacional de Ciencia y Tecnología (CONACYT; to ALA). Partial funding for this work was provided by a grant from CONACyT (2002-C01-39711) to FXG-C. Permits for fieldwork in Mexico and were issued by the Secretaria de Medio Ambiente y Recursos Naturales (SEMARNAT) to FXG-C. We thank the following individuals and institutions for providing tissue loans: R. J. Baker, Natural Sciences Research Laboratory, Texas Tech University; F. A. Cervantes, Universidad Nacional Autónoma de México; M. D. Engstrom, Royal Ontario Museum; R. C. Dower, Angelo State University; and C. Lorenzo, Colección Mastozoológica de El Colegio de la Frontera Sur. Data analyses were performed using the Fulton Supercomputer at Brigham Young University. We also thank the following people for their help with specimen collection, laboratory work, or other critical input: E. Arellano, R. Mercado, D. K. Hardy, and J. L. T. Magadán.
Specimens examined.—For each voucher specimen of Handleyomys we list the museum acronym and catalog number as follows: ASNHC = Angelo State Natural History Collections; BYU = Monte L. Bean Life Science Museum, Brigham Young University; CMC = Colección de Mamíferos del Centro de Investigación en Biodiversidad y Conservación, Universidad Autónoma del Estado de Morelos; CURN = Centro Universitario Regional del Norte de la Universidad Autónoma de Nicaragua; ECOSCM = El Colegio de la Frontera Sur; MZFC = Museo de Zoología, Facultad de Ciencias, Universidad Nacional Autónoma de México; ROM = Royal Ontario Museum; and TCWC = Texas Cooperative Wildlife Collection, Texas A&M University. For sequences from Genbank we list the accession ID. Specimens are listed by taxon, country, collecting location, locality number, museum voucher number, and specimen field number. Abbreviations Cytb, and Fgb-I7 indicate which gene or gene segment was sequenced for each individual.
Handleyomys alfaroi.—ECUADOR: Esmeralda, Comuna San Francisco de Bogotá (37; Cytb = EU579488); MEXICO: Veracruz (Ver), Catemaco, 13.0 km NW (by road) Sontecomapan, Estación Los Tuxtlas-IBUNAM, 150 m (31; CMC 2246 = DSR 8543 [Cytb = KF658401, Fgb-I7 = KF658443], CMC 2247 = DSR 8544 [Cytb = KF658400, Fgb-I7 = KF658444]); NICARAGUA: Matagalpa, Selva Negra, 250 m (27; Cytb = EU579489).
Handleyomys chapmani.—MEXICO: Tamaulipas (Tamps), El Cielo, San José, 1,329 m (1; TCWC 59291 = ICA 36 [Cytb = KF658365, Fgb-I7 = KF658451], TCWC 59294 = ICA 69 [Cytb = KF658373, Fgb-I7 = KF658452], TCWC 59289 = ICA 75 [Cytb = KF658375, Fgb-I7 = KF658450]; San Luis Potosí (SLP), El Naranjo, 3.5 km N 3 km W, Maguey de Oriente (2; CMC 739 = FXG 527 [Cytb = KF658376, Fgb-I7 = KF658422], CMC 740 = FXG 528 [Cytb = KF658356, Fgb-I7 = KF658423], CMC 741 = FXG 529 [Cytb = KF658377); Veracruz (Ver), Zacualpan (3; MZFC 8304 = HBR 069 [Cytb = KF658379, Fgb-I7 = KF658448]); Hidalgo (Hgo), 26.5 km NE (by road) Metepec, 2,210 m (4; CMC 1042 = FXG 804 [Cytb = KF658348], CMC 1043 = FXG 823 [Cytb = KF658353], CMC 1044 = FXG 827 [Cytb = KF658361, Fgb-I7 = KF658431]); Puebla (Pue), Huauchinango, Rancho El Paraíso, 6 km SW Huahuchinango, 2000 m (5; BYU 15801 = EAA 643 [Cytb = KF658362, Fgb-I7 = KF658417], BYU 15802 = EAA 644 [Cytb = KF658354); Puebla (Pue), La Gloria Falls, Apulco River, 10 km N Zacapoaxtla, 1,500 m (6; BYU 15803 = EAA 642 [Cytb = KF658344, Fgb-I7 = KF658418]); Puebla (Pue), 4.7 km NE (by road) Teziutlán, 1,750 m (7; CMC 1049 = FXG834 [Cytb = KF658345, Fgb-I7 = KF658432], CMC 1052 = FXG 837 [Cytb = KF658349], CMC 1054 = FXG 839 [Cytb = KF658346]); Veracruz (Ver), Xico, Matlalapa, 2,070 m (8; CMC 1497 = RMV 50 [Cytb = KF658378], CMC 1495 = RMV48 [Cytb = KF658355, Fgb-I7 = KF658436]); Veracruz (Ver), Xalapa, El Haya, Old road to Coatepec km 25 (Botanic Garden Francisco Javier Clavijero), 1,235 m (9; CMC 1450 = RMV 01[Cytb = KF658343]; Veracruz (Ver), Acajete, Mesa de la Yerba, 3.4 km intersection to Mazatepec (Xalapa-Perote by road), 2,004 m (10; CMC 1353 = FXG 873 [Cytb = KF658366], CMC 1490 = RMV 84 [Cytb = KF658380, Fgb-I7 = KF658435]); Veracruz (Ver), Huatusco, Las Cañadas, 1,340 m (11; CMC 779 = FXG 618 [Cytb = KF658350, Fgb-I7 = KF658426], CMC 780 = FXG 619 [Cytb = KF658360], CMC 782 = FXG 621 [Cytb = KF658351]); Veracruz (Ver), Texhuacán, 1.2 km SE Xochititla, 1,670 m (12; CMC 772 = FXG 578 [Cytb = KF658358, Fgb-I7 = KF658424], CMC 773 = FXG 579 [Cytb = KF658357], CMC 774 = FXG 580 [Cytb = KF658372], CMC 775 = FXG 581 [Cytb = KF658347, Fgb-I7 = KF658425]); Oaxaca (Oax), Puerto de la Soledad, 2,600 m (13; BYU 15303 = EAA 310 [Cytb = KF658364], BYU 15304 = EAA 311 [Cytb = KF658363]); Oaxaca (Oax), Concepción Pápalo, 14.4 km NE (by road) Santa Flor, 2,600 m (14; CMC 1382 = FXG 943 [Cytb = KF658370], CMC 1347 = FXG 944 [Cytb = KF658367], CMC 1352 = FXG 949 [Cytb = KF658381, Fgb-I7 = KF658434], CMC 1389 = FXG 950 [Cytb = KF658369]); Oaxaca (Oax), Ixtlán, 11 km SW (by road) La Esperanza, 2,400 m (15; CMC 113 = DSR 5800 [Cytb = KF658374, Fgb-I7 = KF658419], CMC 115 = DSR 5827 [Cytb = KF658371]); Oaxaca (Oax), Santa María Tlahuitoltepec, Santa María Yacochi, 2,400 m (16; CMC 114 = DSR 5701 [Cytb = KF658368], CMC 117 = DSR 5763 [Cytb = KF658382], CMC 119 = DSR 5765 [Cytb = KF658359]); Oaxaca (Oax), Candelaria Loxicha, 0.7 km E (by road) La Soledad, 1,025 m (17; CMC 943 = FXG 682 [Cytb = KF658395]); Oaxaca (Oax), Miahuatlán, San Miguel Suchixtepec, Río Molino, 2,353 m (18; CMC 925 = FXG 691 [Cytb = KF658388, Fgb-I7 = KF658427], CMC 930 = FXG 737 [Cytb = KF658389], CMC 931 = FXG 738 [Cytb = KF658391, Fgb-I7 = KF658429], CMC 932 = FXG 739 [Cytb = KF658387], CMC 927 = FXG 734 [Cytb = KF658390, Fgb-I7 = KF658428]); Guerrero (Gro), Malinaltepec, 3 km E El Tejocote, 2,620 m (20; CMC 1656 = FXG 1043 [Cytb = KF658392], CMC 1657 = FXG 1044 [Cytb = KF658394, Fgb-I7 = KF658438], CMC 1655 = FXG 1041 [Cytb = KF658393, Fgb-I7 = KF658437]); Guerrero (Gro), Chilpancingo de los Bravos, 6.1 km SW (by road) Omiltemi, 2,480 m (21; CMC 455 = FXG 412 [Cytb = KF658399]); Guerrero (Gro), Leonardo Bravo, 3.4 km (by road) Carrizal, 2,480 m (22; CMC 452 = FXG 462 [Cytb = KF658397, Fgb-I7 = KF658420], BYU 20647 = FXG 463 [Cytb = KF658396], CMC 454 = FXG 464 [Cytb = KF658398, Fgb-I7 = KF658421]); Hidalgo (Hgo), Tlanchinol, 3 km E (by road) Tlanchinol, 1,451 m (23; BYU 15300 = EAA 272 [Cytb = KF658352]).
Handleyomys melanotis.—MEXICO: Oaxaca (Oax), Putla Villa de Guerrero, 5.5 km S (by road) Concepción de Guerrero, 936 m (19; CMC 942 = FXG 789 [Cytb = KF658412, Fgb-I7 = KF658430], CMC 939 = FXG 793 [Cytb = KF658413]; Nayarit (Nay), Peñita de Jaltemba, 1.8 km N of La Peñita de Jaltemba (ASNHC 3418 = ASK1601 [33; Cytb = KF658408]); Michoacán (Mich), Coalcomán, 10.9 km NW (by road) Coalcomán (29; CMC 1806 = DSR 7715 [Cytb = KF658410]); Jalisco (Jal), San Sebastián, 3.4 km W (by road) San Sebastián del Oeste, 1,450 m (30; CMC 1207 = DSR 7414 [Cytb = KF658411, Fgb-I7 = KF658433]); Nayarit (Nay), 8 KM E of San Blas (34; ASNHC 3419 = ASK 1538 [Cytb = KF658409, Fgb-I7 = KF658415]); Colima (Col), Comala, Hacienda San Antonio (36; ASNHC = ASK1957 [Cytb = KF658414, Fgb-I7 = KF658416]).
Handleyomys rostratus.—MEXICO: Veracruz (Ver), Catemaco, 13 km NW (by road) Sontecomapan, Estación Los Tuxtlas, IBUNAM, 150 m (31; CMC 2222 = DSR 8560 [Cytb = KF658407, Fgb-I7 = KF658439]); Chiapas (Chis), Berriozabal, 12 km N (by road) Berriozabal, 1,060 m (32; CMC 2241 = DSR 8464 [Cytb = KF658403, Fgb-I7 = KF658440], CMC 2242 = DSR 8465 [Cytb = KF658406], CMC 2243 = DSR 8466 [Cytb = KF658402], CMC 2244 = DSR 8467 [Cytb = KF658405, Fgb-I7 = KF658441], CMC 2245 = DSR 8468 [Cytb = KF658404, Fgb-I7 = KF658442]); Tamaulipas (Tamps), Rancho Calabazas (near Ciudad Victoria), 3.2 km W Calabazas (35; Cytb = EU579492). EL SALVADOR: Ahuachapán, Ahuachapán, El Imposible (25; Cytb = EU579493). NICARAGUA: Matagalpa, Matagalpa, El Tigre (27; Cytb = EU579491).
Handleyomys saturatior.—MEXICO: Chiapas (Chis), La Trinitaria, Lagos de Montebello (24; ECOSCM 1228 [Cytb = KF658384, Fgb-I7 = KF658446], ECOSCM 1229 [Cytb = KF658385], ECOSCM 1231 [Cytb = KF658383, Fgb-I7 = KF658447]). NICARAGUA: Matagalpa, Selva Negra-Atajo Trail (27; TTU 101644 [Cytb = DQ224410, Fgb-I7 = KF658453]); (28; CURN = JAGE 438 [Cytb = KF658386, Fgb-I7 = KF658445]). EL SALVADOR: Santa Ana, Montecristo National Park (26; ROM 101537 [Cytb = EU579494, Fgb-I7 = KF658449]).
Genealogical sorting index (GSI) and Bayesian phylogenetics and phylogeography (BPP) posterior probability (BpP) for Handleyomys chapmani as currently recognized (Sierra Madre Oriental–Oaxacan Highlands [SMO-OH] and Sierra Madre del Sur [SMS] clades labeled as H. chapmani), and for the SMO-OH and SMS clades labeled as different groups. Values for H. saturatior also are shown as a reference for a recognized and diagnosable species in the group. GSI values correspond to maximum likelihood (ML) and Bayesian inference (BI) topologies from individual genes (GSICytb and GSIFgb-I7), the concatenated data set (GSIConcatenated), and an ensemble from the Cytb and Fgb-I7 trees (GSIT). All GSI values had highly significant P-values (< 0.0004).