Specimens examined and genes sequenced

A total of 116 specimens representing 7 Heteromys anomalus, 3 H. australis, 83 H. desmarestianus, 3 H. gaumeri, 4 H. goldmani, 4 H. nelsoni, 2 H. nubicolens, 6 H. oresterus, 2 Liomys irroratus, and 2 L. salvini was used in this study (Appendix I). Tissue samples were collected from natural populations, or obtained via tissue loans, in Belize, Costa Rica, Guatemala, El Salvador, Honduras, México, Nicaragua, Panama, and Venezuela (Fig. 1) following standard sampling methods approved by the American Society of Mammalogists (Gannon et al. 2007). For this study we used a rapidly evolving mitochondrial gene (cytochrome b [Cytb]) that provided phylogenetic signal for both intra- and interspecific relationships in previous studies of heteromyine rodents (Anderson and Jansa 2006; Rogers and Vance 2005). We also included 2 more slowly evolving nuclear gene segments; engrailed 2 (EN2) exon 3 and myosin heavy polypeptide 6 cardiac muscle alpha (MYH6) exon 35 and intron 35. These latter 2 loci were used to assess deeper nodes in an evaluation of phylogenetic relationships among species of Neotoma (Matocq et al. 2007). In addition to the sequence data generated in this study, Cytb sequences for 1 specimen of H. anomalus and 2 specimens of H. desmarestianus were taken from Rogers and Vance (2005; GenBank accession numbers DQ168468, DQ168466, and DQ168467, respectively), 1 specimen each of H. desmarestianus and H. nubicolens from Anderson and Jansa (2006; DQ450094 and DQ450090, respectively), and 1 specimen of H. gaumeri from Montgelard et al. (2002; AJ389536). L. irroratus and L. salvini were used as outgroup taxa (Hafner et al. 2007) in all analyses. Sequence data for these species were obtained from Rogers and Vance (2005; DQ168477, DQ168483, DQ168491, DQ168494, DQ168497, and DQ168501 for L. irroratus and DQ168540, DQ168542, DQ168545 for L. salvini). Tissue samples of H. oasicus and H. teleus were not available.

DNA extraction, amplification, and sequencing

The entire mitochondrial Cytb gene (1,140 base pairs [bp]), EN2 (189 bp), and MYH6 (252 bp) were amplified and sequenced in both strand directions for this study. Total genomic DNA was extracted from liver, heart, or muscle tissue (either frozen or preserved in 95% ethanol) using the Qiagen DNeasy Tissue Kit (Qiagen Inc., Valencia, California). DNA amplification was accomplished using the polymerase chain reaction (Saiki et al. 1988). Four primers were used to amplify or sequence, or both, the Cytb gene in approximately 800-bp segments: L14724 and H15915 (Irwin et al. 1991), CB3H (Palumbi 1996), and F1 (Whiting et al. 2003). Primers MVZ-16 (Smith and Patton 1993), WDRAT 400F (Edwards et al. 2001), H15149 (Irwin et al. 1991), and Neo700L (Peppers and Bradley 2000) were used as necessary to amplify or sequence, or both, smaller segments. Primers for EN2 and MYH6 were those of Lyons et al. (1997). These authors designed primers for MYH2; however, according to Matocq et al. (2007), the sequences obtained using these primers should be regarded as MYH6.

Parameters for polymerase chain reactions were as follows. For Cytb, 1 cycle of 94°C (3 min) was followed by 36 cycles of 94°C (1 min) denaturing, 46°C annealing (1 min), and 72°C (1 min) extension; the polymerase chain reaction was concluded by 1 cycle of 72°C (3 min). For EN2, 1 cycle of 94°C (10 min) was followed by 32 cycles of 94°C (1 min) denaturing, 57°C annealing (1 min), and 72°C (1 min) extension. For MYH6, 1 cycle of 94°C (10 min) was followed by 32 cycles of 94°C (1 min) denaturing, 62°C annealing (1 min), and 72°C (1 min) extension. Negative (no DNA) controls were run with all amplifications to reveal instances of DNA contamination. Polymerase chain reaction products were visualized on an agarose gel with ethidium bromide. Successfully amplified products were purified with silica gel using the Gene Clean III Kit or by using a Multiscreen PCR 96-Well Filtration System (Millipore Corp., Billerica, Massachusetts). The purified polymerase chain reaction products were cycle-sequenced using the primers described above, and sequenced products purified using a Sephadex protocol or Multiscreen Filter Plates for High Throughput Separations (Millipore Corp.). Light and heavy strand sequences were collected on an ABI 377 automated sequencer (Applied Biosystems, Foster City, California) and then edited and compiled using Sequencher versions 3.1.1 and 4.1.2 (Gene Codes Corp., Ann Arbor, Michigan). Alignments for Cytb and EN2 were unambiguous (no insertions–deletions [indels]) and were performed using Se-Al version 2.0 (Rambaut 2002). MYH6 was aligned using Clustal_X (Thompson et al. 1997). Sequences were submitted to GenBank (GenBank accession numbers GU646919–GU647038 for Cytb, GU731466–GU731514 for EN2, and GU657039–GU647048 and GU731427–GU731465 for MYH6; Appendix I).

Phylogenetic analyses of the Cytb data set

From the total Cytb data set of 116 individuals COLLAPSE version 1.2 (Posada 2004) identified 77 unique haplotypes. Phylogenies were estimated for these 77 haplotypes using maximum parsimony (MP) as implemented in PAUP* 4.0b10 (Swofford 2003), maximum likelihood (ML) as implemented in PhyML (Guindon and Gascuel 2003), and Bayesian inference (BI) with Markov chain Monte Carlo sampling for mixed models using MrBayes 3.0b4 (Ronquist and Huelsenbeck 2003). Model selection (both linked and unlinked) for BI was performed for each codon position using the Akaike information criterion (AIC—Akaike 1974) as implemented in MrModeltest version 3.7 (Nylander 2004) and as suggested by Posada and Buckley (2004). The best-fit model of evolution for the 1st and 3rd positions was TrN+I+Γ and HKY+Γ for the 2nd position. Bayesian posterior probabilities (pP) were determined by running 20 chains for 20 million Metropolis coupled Markov chain Monte Carlo generations using the default priors on model parameters. For all analyses 1,000 trees were sampled from the posterior probability distribution (1 every 1,000 generations). A majority-rule consensus tree was produced after discarding the burn-in determined by Tracer version 1.4 (Rambaut and Drummond 2003; available from  http://tree.bio.ed.ac.uk/software/tracer/). Each BI analysis was conducted twice, starting from a different, randomly chosen tree. Posterior probabilities for individual clades obtained from independent analyses were compared for congruence (Huelsenbeck and Imennov 2002; Huelsenbeck et al. 2002; Nylander 2004). Percent sequence divergence was estimated for this data set and the combined data set (see below), using both MP and ML as implemented in PAUP* 4.0b10 (Swofford 2003). MP analyses were conducted with equal character weighting with tree-bisection-reconnection branch swapping.

Under the ML criterion the model of evolution most appropriate for the Cytb gene was HKY+I+Γ (Hasegawa et al. 1985), as selected using AIC as implemented in Modeltest version 3.6 (Posada and Crandall 1998). The base frequencies were A  =  0.3204, C  =  0.2897, G  =  0.0816, and T  =  0.3083; transversion rates were (A–C) 1.0000, (A–G) 15.9820, (A–T) 1.6077, (C–G) 1.6077, and (C–T) 15.982. The proportion of invariable sites (I) was 0.5540, and the gamma distribution shape parameter (Γ) was 1.2844.

Phylogenetic analyses of the EN2 and MYH2 data sets

From the MP Cytb haplotype tree a subset of 46 ingroup terminals was selected and sequenced for the EN2 and MYH6 gene segments following the hierarchical reduced sampling protocol of Morando et al. (2003). The evolutionary models selected by Modeltest and MrModeltest were HKY+I+Γ for EN2 and GTR+Γ for MYH6. These models were used in both the ML and BI analyses.

Phylogenetic analysis of combined data set

In all but 4 instances the same specimen was sequenced for all 3 genes for the 46 ingroup terminals. In these 4 instances another individual from the same Cytb haploclade was used as a substitute. With the exception of a single individual from Panama (H. desmarestianus complex clade VII, locality 45; see Appendix I), for which we could not obtain either EN2 or MYH6 sequence, all Cytb haploclades were represented in all combined data analyses.

Following Wiens (1998) and Liu and Miyamoto (1999), we analyzed each gene individually and then with all sequences concatenated into a single analysis to explore the extent of congruence and conflict among of data sets more fully. Because of low numbers of parsimony-informative characters in the EN2 and MYH6 gene segments (González 2005) and the computational time required, MP analyses for these gene segments were not performed. For the MP analysis of combined data, parameters were the same as for the Cytb haplotype data set. We found 2 gaps (each a single base pair in length) resulting from 2 indel events in MYH6. These 2 gaps each were treated as a 5th state. For the BI analysis Cytb partitions or models, or both, and models of evolution for EN2 and MYH6 remained the same.

Under the ML criterion the model of evolution most appropriate for the combined data was (HKY+I+Γ). The base frequencies were A  =  0.2883, C  =  0.2807, G  =  0.1731, and T  =  0.2579; transversion rates were (A–C) 2.3510, (A–G) 10.0432, (A–T) 3.8820, (C–G) 0.9356, and (C–T) 32.5061. The proportion of invariable sites (I) was 0.6124, and the gamma distribution shape parameter (Γ) was 1.2170.

Nodal support

For MP trees branch support (BS) for nodes was estimated using nonparametric bootstrapping (Felsenstein 1985) with 10,000 pseudoreplicates of heuristic searches using tree-bisection-reconnection branch swapping and 100 random sequence additions. Additionally, partitioned Bremer supports (Baker and DeSalle 1997; Baker et al. 1998; Bremer 1988) were calculated using TreeRot, version 2c (Sorenson 1999) to evaluate support for nodes for each gene or gene segment. Negative Bremer support values indicate conflict among gene trees at that particular node. For Bayesian analyses trees not discarded as burn-in were used to construct a 50% majority rule consensus tree, with values representing posterior probabilities. Posterior probabilities (pP) ≥0.95 were considered evidence of significant support for a clade (Huelsenbeck and Ronquist 2001). BS support values for ML, reported as bootstrap proportions (BP), were estimated using nonparametric bootstrapping with 1,000 pseudoreplicates in PhyML (Guindon and Gascuel 2003).

Hypothesis testing

Alternative phylogenetic hypotheses were evaluated for statistical significance using the 1-tailed Shimodaira–Hasegawa test (Shimodaira and Hasegawa 1999) with restricted ML, a method of estimating components of variance as implemented in PAUP* 4.0 (Swofford 2003). Ten thousand bootstrap replicates were performed using the Shimodaira–Hasegawa topology test by resampling the partial likelihoods for each site (RELL model).

Delimiting species boundaries

The discovery of species is generally recognized as 1 of the 2 stated empirical goals of systematic biology (Wiley and Mayden 2000), and the number of species definitions or concepts is numerous (Baker and Bradley 2006; de Queiroz 1998; Wheeler 1999; Wiley and Mayden 2000). Although these definitions typically allude to the same general idea of what a species is, in many instances workers propose different properties as being the essential determinant of a species. Unfortunately, the majority of species definitions provide no framework with which species boundaries actually could be tested. Recently, this area of systematics has received much needed attention (Sites and Marshall 2003, 2004; Wiens 2007). Although the methods may not agree, progress continues toward viewing species limits as testable hypotheses rather than vague assemblages of characters or de facto thresholds of genetic distance (Bradley and Baker 2001). The DNA tree-based approach of species delimitation of Wiens and Penkrot (Wiens and Penkrot 2002) was selected a priori to delimit species in a hypothesis-testing framework. This method uses a dichotomous key approach in making decisions at the species level and takes advantage of the more rapid coalescing time of mtDNA. A focal species is selected and examined for exclusivity ( = monophyly), meaning that haplotypes representing the focal species do not interdigitate with haplotypes representing other entities. If other clusters of haplotypes are recovered as strongly supported basal clades, which are exclusive by locality, those terminals also may be considered a species-level clade. This method can identify cryptic species, particularly when well-supported clades are concordant with geography and may have been overlooked based on morphological comparisons (Wiens and Penkrot 2002).

We selected H. desmarestianus as our focal species, because examination of karyotypic (Rogers 1989) and allozyme (Anderson et al. 2006; Rogers 1990) data suggested a strong possibility of multiple species. The method of Wiens and Penkrot (Wiens and Penkrot 2002) was used to delimit species at the deepest level of divergence in the tree, and we followed the protocol of Morando et al. (2003) in defining basal clades as those represented by at least 2 haplotypes and separated beyond the TCS (Clement et al. 2000) network connections (of 95% confidence; 14 steps for this data set—González 2005) among Cytb haplotypes. After using the approach of Wiens and Penkrot (Wiens and Penkrot 2002), we incorporated data from previous studies of Heteromys (e.g., morphology and karyology) in our determination of species limits. Uncorrected “p” distances (p-distance) for Cytb were generated using PAUP* to allow for comparison of relative sequence divergence among clades but were not used in assessing species-level status.

Phylogenetic analysis of individual genes

Analyses of the Cytb haplotypes data set (1,140 bp) based on MP, ML, and BI optimality criteria converged on essentially the same tree topology, but with less resolution using the MP optimality criterion (results available on request). MP analysis of the Cytb haplotypes data set resulted in a single most-parsimonious tree (length  =  1,944; consistency index  =  0.342; retention index  =  0.818, tree not shown). The data matrix contained 461 variable characters, of which 426 were parsimony informative (González 2005). For the ML tree topology, except for clade II, nodal support was high for all species-level taxa and haploclades (Fig. 2). Of the 77 Cytb haplotypes evaluated in these analyses, all were recovered as exclusive by locality with a few exceptions. A single individual of H. anomalus from Miranda, Venezuela (locality 2 in Fig. 1), grouped with individuals from Sucre (locality 1). Several haplotypes representing H. desmarestianus also were found in more than 1 locality; these instances were confined to several locations in the Guatemalan Peten and the Yucatan Peninsula of Mexico (samples 19 and 21; samples 13–15 and 20).

Figure 1

Map of Mexico and Central America with an inset of northern South America illustrating the geographic location of samples included in this study. Numbers for localities are the same as those used in Figs. 2 and 3 and Appendix I. Samples of Heteromys are represented by • and Liomys by ▴.

Figure 2

Maximum-likelihood phylogenetic hypothesis (log-likelihood score −1,070.58) for Heteromys based on Cytb sequence data using the HKY+I+Γ model of evolution and rooted with Liomys salvini and L. irroratus. Numbers to the right of terminals indicate sampling localities (Appendix I), followed by specimen voucher numbers and abbreviations designating countries in which the samples were collected as follows: BL  =  Belize; CR  =  Costa Rica; ES  =  El Salvador; GU  =  Guatemala; HN  =  Honduras; MX  =  Mexico; NI  =  Nicaragua; PN  =  Panama; VZ  =  Venezuela. Nodal support is represented by maximum-likelihood (ML) bootstrap replicates (1,000 iterations) and bootstrap proportions (BP) above nodes, maximum-parsimony (MP) nonparametric bootstrap pseudoreplicates and values (branch support [BS] ≥70% shown with an asterisk and derived from 10,000 pseudoreplicates) below nodes, and solid circles at nodes indicating Bayesian posterior probabilities (pP) ≥0.95%. Bootstrap support for terminal nodes (samples from a single locality) is not shown. Roman numerals refer to clades currently regarded as H. desmarestianus and as detailed in text. Letters indicate candidate species as discussed in text.

Heteromys anomalus, H. australis, H. gaumeri, and H. nelsoni were recovered as highly supported monophyletic clades (pP ≥ 0.95, BP  =  100, BS ≥ 70%). The remaining specimens grouped together in what we refer to as the H. desmarestianus complex, composed of 11 clades (Fig. 2), 3 of which currently are recognized as H. goldmani, H. nubicolens, and H. oresterus. Within this complex both ML and BI analyses recovered relationships among clades I, VII, and VIII as (VII–VIII) (I); however, this grouping is not well supported and was not present in the MP tree. Relationships among the remaining clades and substructure within clades were largely concordant across optimality criteria.

Bayesian and ML analyses of EN2 and MYH6 resulted in less-resolved trees in comparison with our analysis of Cytb (González 2005). Analyses of EN2 sequences yielded the least-resolved trees, which was not surprising because of few informative characters (13 of 15 variable characters). However, BI and ML analyses of the EN2 data set recovered a clade consisting of a H. desmarestianus complex (clades II, III, IV, and V in Fig. 2) together with H. goldmani, H. nubicolens, and H. oresterus. Additionally, clade I was recovered as monophyletic, as was H. gaumeri. EN2 analyses did not recover the genus Heteromys as monophyletic relative to Liomys, although support for this nonmonophyly was not significant. Bayesian analysis of the MYH6 gene segment resulted in a more-resolved consensus tree than did EN2 (González 2005). The MYH6 gene segment yielded 44 variable characters of which 30 were parsimony informative. In addition, 2 single base pair indels were identified within the MYH6 data set. The 1st defined a clade consisting of Liomys + Heteromys, and the 2nd indel was a synapomorphy for the H. desmarestianus species group sensu Goldman (1911) and Hall (1981), which includes the H. desmarestianus complex (including H. goldmani, H. nubicolens, and H. oresterus). In both BI and ML analyses of the MYH6 gene Heteromys was recovered as monophyletic, as were H. nelsoni and the H. desmarestianus complex clades IV, VII, and H. oresterus (González 2005).

Phylogenetic analysis of combined data set

Using the qualitative method of evaluating congruence among data sets as described by Wiens (1998), Cytb, MYH6, and EN2 tree topologies were compared and instances of incongruence were noted. MYH6 and EN2 topologies were not well resolved, and the majority of nodes in the trees were comprised of polytomies that did not disagree with nodes resolved in the Cytb topology. The only incongruence was in the EN2 topology in which a clade composed of H. australis, H. anomalus, and Liomys irroratus was recovered (pP  =  0.78—data available on request).

Trees estimated from concatenating Cytb, MYH6, and EN2 sequence data converged on essentially the same tree topology regardless of the optimality criterion used; differences were the result of less resolution. The ML tree (lnL  =  −9,765.58; Fig. 3) recovered the same major relationships as the Cytb tree (Fig. 2), with the exception that groups I and VII were not associated with each other. This was due, at least in part, to “clade” VIII not being represented in the combined analysis. Bayesian pP were significant, and ML bootstrap values were high (BP ≥ 86) for all species-level taxa and haploclades. Partitioned Bremer support indicated that most of the support for the topology was from the Cytb data set; however, 9 of 15 interior nodes were supported by 2 or more genes as indicated by positive Bremer support values, including monophyly of the genus Heteromys and monophyly of the H. desmarestianus group. H. anomalus, H. australis, and H. gaumeri each were recovered as monophyletic groups with high support in the combined analysis (pP ≥ 0.95, BP  =  ≥ 99). BI and ML analyses supported a clade composed of H. anomalus and H. australis, but Bremer support for this association was not strong. H. nelsoni consistently was recovered as a basal clade relative to the other species in the genus. The remaining clades in the topology comprise the H. desmarestianus complex, which includes H. goldmani, H. nubicolens, and H. oresterus. This latter clade also was supported by a single base pair deletion in MYH2. We view the ML topology with Bayesian analysis of the combined data set (Fig. 3) to be our best working hypothesis and use this topology as our tree of reference.

Figure 3

Maximum-likelihood phylogenetic hypothesis (log-likelihood score −9,765.58) for the genus Heteromys based on combined Cytb, MYH6, and EN2 sequence data using the HKY+I+Γ model of evolution and rooted with Liomys salvini and L. irroratus as outgroup taxa. Numbers to the right of terminals indicate sampling localities (Appendix I), followed by specimen voucher numbers and country abbreviations as in Fig. 2. Closed bars at nodes indicate unique MYH6 indels. Nodal support from maximum-likelihood and Bayesian posterior probabilities is as in Fig. 2. Partitioned Bremer support values (Cytb/MYH6/EN2) are provided below nodes. Bootstrap support for terminal nodes is not shown. Roman numerals and letters are as defined in Fig. 2.

Inferred species boundaries

Heteromys anomalus, H. australis, H. gaumeri, and H. nelsoni were recovered as strongly supported clades in all analyses, thus supporting their species-level status. Following the dichotomous key presented in Wiens and Penkrot (2002), we determined that our focal species, H. desmarestianus, was nonexclusive (i.e., paraphyletic) with respect to 1 or more distinct, exclusive species (H. goldmani, H. nubicolens, and H. oresterus). Based on criteria of the method of Wiens and Penkrot (Wiens and Penkrot 2002), the H. desmarestianus complex has 11 distinct haplotype groups exclusive by locality (clades I–VIII and H. goldmani, H. oresterus, and H. nubicolens; Fig. 2).

Hypothesis testing

Using the Shimodaira–Hasegawa test (Shimodaira and Hasegawa 1999), we tested 8 a priori hypotheses derived from the literature (Table 2) for the Cytb and the combined data sets. All constraint tests resulted in significantly less-likely tree topologies with the exception of monophyly of the H. anomalus species group for both Cytb and the combined data.

Table 2

Topological tests of a priori hypotheses using Shimodaira–Hasegawa tests (Shimodaira and Hasegawa 1999). Log-likelihood scores for the cytochrome-b and combined data trees were −1,070.58 and −9,765.58, respectively. Probabilities ≤0.05 indicate statistical rejection of the respective hypothesis.

Relationships among species of Heteromys

We tested the relationships developed by Rogers (1990) and Anderson et al. (2006) against the tree topologies depicted in Figs. 2 and 3. The relationships among taxa recovered from Cytb or the combined gene analyses were significantly more likely than the topology based on allozymes (Rogers 1990) or on the alternative topology based on a combination of allozymes and morphology (Anderson et al. 2006). Aside from the topological differences between our phylogenetic hypothesis and those developed earlier, the differences may be due, at least in part, to earlier trees being less resolved because both the allozyme and morphological data sets were composed of fewer characters.

Heteromys nelsoni consistently was recovered as a basal clade relative to the remaining taxa in the genus (Figs. 2 and 3). In his phenetic analysis of morphological and allozyme characters, Rogers (1986, 1990) also recovered H. nelsoni as distinct and basal to all other Heteromys. No support was found for the subgenus Xylomys consisting of H. oresterus and H. nelsoni (Anderson et al. 2006; Hall 1981). Moreover, forcing a sister-group relationship between H. nelsoni and H. oresterus resulted in a significantly less-likely tree (Table 2). Based on the results of Rogers (1986, 1990), the present study and the unique karyotype possessed by H. nelsoni (2n  =  42—Rogers 1989), maintaining the subgenus Xylomys, composed only of H. nelsoni, is supported (test 1; Table 2).

Table 3

Summary of taxonomic recommendations for members of the genus Heteromys included in this study and evidence for species-level status. Karyotypic data are summarized by Anderson et al. (2006), Patton and Rogers (1993), and Rogers (1989); allozyme data are from Rogers (1986, 1990); and morphological data are summarized by Rogers and Schmidly (1982), Anderson and Jarrín-V. (2002), Anderson and Timm (2006), and Anderson et al. (2006). Habitat preference is derived from museum collection records and Anderson (1999, 2003).

The H. anomalus species group is represented here by H. anomalus and H. australis and represents a 2nd, basal clade. This arrangement is supported by combined cladistic analysis of morphological and allozyme data by Anderson et al. (2006) but is discordant with allozyme results of Rogers (1990), which revealed no affinity or similarity between H. anomalus and H. australis. Forcing nonmonophyly of the H. anomalus species group did not result in significantly less-likely trees (test 4; Table 2; P  =  0.082). H. australis is restricted in its distribution, known only from humid evergreen forests in eastern Panama, western and Andean Colombia, northwestern Venezuela, and northwestern Ecuador (Anderson 1999, 2003), whereas H. anomalus occurs in deciduous and evergreen tropical forests in Colombia, Venezuela, Trinidad, and Tobago (Anderson 2003). Although the test for nonmonophyly was not significant, evidence from this study does not refute the sister-group relationship between H. anomalus and H. australis. Moreover, other South American taxa (H. oasicus, H. teleus, and H. sp.) were not included in our analysis. According to Anderson et al. (2006), H. oasicus and H. teleus share morphological features with H. anomalus and H. australis. Therefore, we follow Anderson et al. (2006) by placing both taxa in the H. anomalus species group (Table 3). In addition, Anderson and Gutiérrez (2009) described H. catopterius from Venezuela and indicated that this taxon corresponded to part of the samples considered H. anomalus by Anderson et al. (2006). Therefore, we tentatively place this new species in the H. anomalus species group.

Individuals of H. anomalus included in our analyses did not sort by locality. Specifically, a specimen from Miranda, Venezuela (CM 78170), formed a clade with specimens from Sucre. These results were concordant with those of Rogers (1986), who found 3 fixed allozyme differences between mice from Sucre and most from Miranda; however, CM 78170 from Miranda shared the same alleles with individuals from Sucre. The p-distances for Cytb between these 2 clades range from 7.8% to 8.0%. This p-distance is comparable to the 7.7% sequence divergence between the sister taxa Liomys adspersus (Panamanian spiny pocket mouse) and L. salvini (Rogers and Vance 2005). Additional sampling within Venezuela is needed to determine if the 2nd clade we recovered represents H. catopterius or another species-level entity within H. anomalus proper.

Heteromys gaumeri was recovered as the sister clade to the H. desmarestianus complex in all analyses, with strong nodal support. However, H. gaumeri showed no well-supported kinship with H. desmarestianus in allozyme or morphological analyses (Anderson et al. 2006; Rogers 1986, 1990) and has a unique karyotype (2n  =  56, FN  =  76). Based on allozyme, morphological, and chromosomal differences together with results of our study, we support the recommendation of Engstrom et al. (1987) that H. gaumeri be removed from the desmarestianus species group and placed in a group of its own.

Species-level phylogenetics in the H. desmarestianus group

As currently recognized, H. desmarestianus is not monophyletic. Nodal support for this complex (8 clades previously referable to H. desmarestianus together with H. goldmani, H. nubicolens, and H. oresterus) is significant in terms of Bayesian posterior probabilities but is not particularly strong based on MP or ML bootstrap values. This complex also is defined by a single base pair deletion in MYH6. Forcing nonmonophyly of the H. desmarestianus group resulted in significantly less-likely trees for both Cytb and the combined data sets (test 5; Table 2).

The tree topology of the H. desmarestianus species group has a high level of geographic structure (Fig. 4). Clades I, II, III, and H. goldmani are located in southern Mexico (states of Veracruz, Oaxaca, and Chiapas, respectively). Clade IV is composed of specimens restricted to Central America (Mexican states of Campeche and Quintana Roo, Belize, Guatemala, El Salvador, and Honduras); clade V is found in Costa Rica and Nicaragua; H. oresterus, H. nubicolens, and clade VI are located in Costa Rica; clade VII occurs in eastern Costa Rica and western Panama; and clade VIII is made up exclusively of specimens from the Panamanian Darien.

Figure 4

Map of Mexico and Central America illustrating the geographic location of species, candidate species, and major clades within the Heteromys desmarestianus complex. Roman numerals and letters are those used in Figs. 2 and 3 and the text.

Key to sorting out the nonexclusivity of H. desmarestianus is the question of whether H. goldmani, H. nubicolens, and H. oresterus merit species-level status. Otherwise, the H. desmarestianus complex could be viewed as a widespread, monophyletic entity with high levels of genetic and karyotypic diversity. Earlier studies all have documented subdivisions within H. desmarestianus (Anderson et al. 2006; Rogers 1986, 1989, 1990), and with the range of H. desmarestianus spanning nearly 2,000 km from Mexico to eastern Panama and northern South America (Rogers 1986; Williams et al. 1993), geographic extent alone would argue for the possibility of cryptic species. Inasmuch as H. goldmani is nested more internally than either H. nubicolens or H. oresterus (Figs. 2 and 3), we address this taxon first.

Rogers and Schmidly (1982) noted that compared to H. d. desmarestianus, H. goldmani was darker and lacked a pronounced sprinkling of ochraceous hairs, with several larger cranial features and a smaller, more rounded baculum. Rogers (1986) included additional samples of H. desmarestianus and found that although H. goldmani averaged larger than H. desmarestianus in many cranial features, it fell within the range of morphological variation present in H. desmarestianus from southern Mexico. Rogers (1990) determined that H. goldmani was not significantly different morphologically from the adjacent populations of H. desmarestianus, but he also documented 2 fixed allozyme differences relative to H. desmarestianus. In all of our analyses H. goldmani was recovered as a strongly supported, monophyletic clade (Figs. 2 and 3) separate from nearby H. desmarestianus in Mexico (Fig. 4). H. goldmani is genetically divergent from H. desmarestianus (p-distances ranging from 8.0% to 9.2%) and possesses a unique karyotype (2n  =  60, FN  =  78—Rogers 1989). Constraint tests for the monophyly of H. goldmani relative to adjacent Mexican populations of H. desmarestianus result in significantly less likely trees for a nonmonophyletic H. goldmani (test 6; Table 2). According to the flow chart for species delimitation (Wiens and Penkrot 2002), Goldman's spiny pocket mouse fits case “1c.” Therefore, we regard H. goldmani as a species-level taxon and recommend that it remains in the H. desmarestianus species group within the subgenus Heteromys.

Heteromys oresterus (San José and Cartago provinces, Cordillera Central, Costa Rica) consistently forms a sister-group relationship with individuals of H. nubicolens from the nearby provinces of Guanacaste and Puntarenas (Cordillera de Guanacaste and Tilarán; Figs. 2 and 3, localities 40 and 41). This relationship was not recovered by Rogers (1986, 1990), who included a sample of H. nubicolens from Monte Verde, Guanacaste Province (then regarded as H. desmarestianus), in his allozyme analysis. Likewise, a cladistic analysis of allozyme and morphological data for members of the subfamily Heteromyinae by Anderson et al. (2006) did not recover H. nubicolens and H. oresterus as sister taxa. However, in contrast to support values in those analyses, nodal support for the sister-group relationship between these 2 taxa is strong across all of our data sets and analyses. H. oresterus (FN  =  78) differs from H. nubicolens (FN  =  86) karyotypically and allozymically (7 fixed differences—Rogers 1990), and Cytb p-distances between these 2 taxa range from 6.6% to 6.9%. Qualitative morphological differences also exist (Anderson and Timm 2006). Finally, forcing nonmonophyly of samples regarded as either H. oresterus or H. nubicolens results in significantly less likely trees (tests 7 and 8; Table 2). Both taxa conform to case “1e” of the method of Wiens and Penkrot (Wiens and Penkrot 2002). Therefore, we recognize these 2 taxa as valid sister species belonging to the H. desmarestianus species group (Table 3).

The clade composed of H. oresterus and H. nubicolens is divergent genetically from H. desmarestianus, clade VI, which is comprised of spiny pocket mice from the Caribbean coastal province of Limón (Fig. 3) and 1 individual from the inland province of Alajuela (Fig. 2). Both H. oresterus and H. nubicolens occur in cloud forests, whereas individuals representing clade VI inhabit low- or midelevation forests (Table 3). Although the nondifferentially stained karyotypes of both Heteromys from coastal Limón and H. nubicolens are FN  =  86 (cytotype F in Rogers [1989]), the Cytb p-distances between the H. oresterus–H. nubicolens clade and clade VI range from 10.0% to 10.3%. Likewise, spiny pocket mice from Limón differ by 3 and 7 allozyme fixed differences from individuals representing H. nubicolens and H. oresterus, respectively (Rogers 1990). These results fit case “1e” of the approach of Wiens and Penkrot (Wiens and Penkrot 2002). Therefore, we consider clade VI to represent candidate species A (Figs. 2 and 3).

Clade I forms the sister group to clades II–VI together with H. goldmani, H. nubicolens, and H. oresterus (Fig. 3). Spiny pocket mice representing this clade are known from low-elevation localities in Oaxaca and Veracruz, Mexico (Fig. 4), and possess an FN  =  82 or FN  =  86 karyotype (Rogers 1989). This clade was recognized by Rogers (1990), who documented 2 allozyme fixed differences distinguishing these lowland populations from other samples of H. desmarestianus in Mexico and Honduras. Clade I also was recovered by Anderson et al. (2006) in their analysis of allozyme and morphological data for the subfamily Heteromyinae. Clade I is highly divergent genetically (Cytb p-distances range from 11.6% to 15.8%) from all other Heteromys, including other samples of H. desmarestianus. This pattern is consistent with case “1c” of the method of Wiens and Penkrot (Wiens and Penkrot 2002). We refer to this entity as candidate species B (Figs. 2 and 3).

Based on the combined data analyses, clade VII is basal to the rest of the H. desmarestianus complex. The Costa Rican and Panamanian samples forming this clade separated into 2 haplotype networks (González 2005), and the Cytb p-distance between these 2 subclades was 8.7%. This clade is well supported in all analyses (pP  =  1.0, BP  =  100). Compared with other Heteromys in our analyses, Cytb p-distances for clade VII range from 12.8% to 17.0%. The Costa Rican individuals included in this study were identified as Heteromys sp. by Mascarello and Rogers (1988) and H. desmarestianus planifrons by Anderson et al. (2006). Mascarello and Rogers (1988) documented standard and differentially stained karyotypic differences (Table 3) between these samples and both H. desmarestianus, from near Tilarán, Province Guanacaste, Costa Rica (clade V), and H. oresterus. Compared with other Heteromys examined by Rogers (1990), the Costa Rican samples of clade VII were distinguished by 2 fixed allozyme differences. This taxon fits case “1e” of the method of Wiens and Penkrot (Wiens and Penkrot 2002). We consider this clade to be candidate species C (Figs. 2 and 3).

Clade VIII is composed of a single individual from Cana, Panama (Fig. 4). Although collected at a locality proximate to H. australis, this individual is highly divergent from all other Heteromys (Cytb p-distances ranged from 10.9% to 15.8%), was identified by Anderson et al. (2006) as H. d. crassirostris, and differs from H. australis morphologically (Anderson 1999). Unfortunately, this animal was not included in the combined analyses because of our inability to obtain nuclear sequence data. However, this individual also possesses 2 fixed allozymic differences (Rogers 1990) compared with samples of H. anomalus, H. australis, or H. desmarestianus. Although the method of Wiens and Penkrot (Wiens and Penkrot 2002) could not be applied to this entity because it is represented by a single individual, we consider this clade to be candidate species D (Fig. 2).

Within the “interior” of the H. desmarestianus complex (see Figs. 2–Figure 34), there are 2 reciprocally monophyletic groups that each exhibit substructure. The 1st (clades II, III, and IV) is found in Mexico, Belize, El Salvador, and Guatemala with low FN values (67, 68, or 72). Members of this clade inhabit low- to moderate-elevation tropical forests. The 2nd group (clade V) is mostly Costa Rican and includes 2 localities with an FN  =  86 karyotype together with a single individual from Nicaragua. This second group also occupies both low- and moderate-elevation habitats. These low and high FN groups formed independent haplotype networks at the 95% confidence level (González 2005), but the average Cytb p-distance between clades II–IV and clade V is moderate (4.9–6.3%). Because our study included only 2 samples of H. desmarestianus from the region between these 2 groups (Honduras, clade IV and Nicaragua, clade V), it is possible that more sampling within Central America could provide the haplotypes necessary to connect these 2 clades. However, given the level of karyotypic variation, moderate genetic differentiation, and diversity in habitat preference, we predict that multiple species-level entities exist within these groups.

Recognizing clades for candidate species status has been aided greatly by the advance of DNA sequence data, because it is now possible to find divisions within species that have been overlooked based on gross morphological comparisons. However, this information must be used with caution because DNA sequences sometimes can make it easy to distinguish populations so that some workers feel justified in declaring unique haplotypes a species (Agapow et al. 2004). In addition, incomplete lineage sorting and reliance on a single marker (or linked markers) can lead to spurious phylogenetic reconstructions (Moore 1995). Therefore, we have taken a conservative approach in defining candidate species by following a protocol of tree-based species delimitation and seeking concordance with evidence provided by earlier studies of Heteromys. Although use of genetic data is a powerful tool, it is essential to ensure that species designations retain both meaningful and practical value, and this requires incorporating as much supporting data as possible. Unfortunately, the majority of characters helpful in delimiting some species of Heteromys involve genetic features rather than discrete and easily discernable external morphological characters that are most useful to field workers.

Our recommendations for Heteromys taxonomy are summarized below and in Table 3. We maintain the subgenera Heteromys and Xylomys, with the latter remaining monotypic (H. nelsoni). We divide the subgenus Heteromys into 3 species groups: the H. anomalus group (H. anomalus and H. australis), H. desmarestianus group (H. desmarestianus, H. goldmani, H. nubicolens, H. oresterus, and candidate species A, B, C, and D), and the H. gaumeri group (H. gaumeri). We tentatively assign H. teleus and H. oasicus to the H. anomalus species group following Anderson et al. (2006) and suggest that a new species of Heteromys described by Anderson and Gutiérrez (2009) also should tentatively be placed in the H. anomalus group.

Valid names for candidate species cannot be suggested at this time, inasmuch as many names are available either in synonymy or as subspecies of H. desmarestianus, and determining the valid name (or proposing a new one) will require additional investigation. Based on their geographic proximity to the localities representing candidate species A, several subspecific names might be applied including H. d. planifrons (Goldman 1937), H. d. underwoodi (Goodwin 1943), and H. d. zonalis (Goldman 1912). Specimens representing candidate species B include individuals collected from the type locality (locality 7 in this study) of H. temporalis (Goldman 1911), a name in synonymy with H. desmarestianus. However, other names with priority might properly be applied to this taxon, including H. longicaudatus, with a type locality of “Mexico” (Gray 1868) and defined by Goldman (1911) to match most nearly specimens of Heteromys from Montecristo [Emiliano Zapata], Tabasco, and H. lepturus (Merriam 1902) described from near Guichicovi, Oaxaca. Genetic data from topotypes would shed light on which, if any, of the available names for this candidate species should be applied. Candidate species C is known from the Pacific lowlands in southern Costa Rica and northern Panama (Fig. 4). Here, the 2 most likely available names are H. d. chiriquensis (Enders 1938) and H. d. repens (Bangs 1902). Candidate species D is represented by only a single specimen from the Panamanian Darien. The most probable name to be applied to this entity would be H. d. crassirostris (Goldman 1912). Clearly, the relatively dense geographic sampling of museum specimens in Mexico and Central America compared with South America will facilitate the morphological revisionary studies necessary to assign names to these candidate taxa.