Introduction

In North and Central America, mice in the genus Peromyscus are among the most frequently encountered species and are considered one of the small mammals better studied in many biological subdisciplines (^{Bradley et al., 2004}, ²⁰⁰⁷; ^{Sullivan et al., 2017}). However, fundamental aspects of their biology, such as their species-level diversity, are still unclear (^{Bradley et al., 2007}; ^{Miller & Engstrom, 2008}; ^{Platt et al., 2015}; ^{Sullivan et al., 2017}). For example, the number of recognized Peromyscus species has shifted recently from 66 (^{Pardiñas et al., 2017}) to 76 (^{Álvarez-Castañeda et al., 2019}; ^{Bradley et al., 2019}; ^{Greenbaum et al., 2019}; ^{León-Tapia et al., 2020}; ^{López-González et al., 2019}), and further evidence suggests additional undescribed diversity (^{Ávila-Valle et al., 2012}; ^{Bradley et al., 2017}; ^{Castañeda-Rico et al., 2014}; ^{Durish et al., 2004}; ^{Pérez-Consuegra & Vázquez-Domínguez, 2017}). Moreover, Peromyscus is a paraphyletic taxon because other phenotypically distinct genera (Habromys, Megadontomys, Neotomodon, Osgoodomys, and Podomys) are phylogenetically nested within Peromyscus (^{Bradley et al., 2007}; ^{Miller & Engstrom, 2008}; ^{Platt et al., 2015}; ^{Sullivan et al., 2017}). As such, the current taxonomic arrangement of Peromyscus probably provides only a rough approximation of their species diversity and systematics.

Morphometrical, morphological, karyological, and biochemical data were amply used for nearly a century to understand phylogenetic relationships of Peromyscus mice (^{Carleton, 1980}, ¹⁹⁸⁹; ^{Hsu & Arrighi, 1966}; ^{Osgood, 1909}; ^{Stangl & Baker, 1984}), and these early findings were subsequently refined with the use of DNA sequencing (^{Bradley et al., 2007}; ^{Miller & Engstrom, 2008}; ^{Platt et al., 2015}; ^{Sullivan et al., 2017}). Most of these more recent investigations sequenced the mitochondrial gene cytochrome b (cyt-b), which has proven beneficial to clarify phylogenetic relationships between species in Peromyscus and related genera, such as Baiomys, Habromys, Megadontomys, Nelsonia, Neotoma, and Reithrodontomys (^{Arellano et al., 2005}; ^{Bradley et al., 2017}, ²⁰¹⁹; ^{Hernández-Canchola et al., 2021}; ^{Hernández-Canchola & León-Paniagua, 2021}; ^{León-Tapia, 2013}; ^{Rogers et al., 2007}; ^{Vallejo & González-Cózatl, 2012}).

Within Peromyscus, species groups have provided valuable taxonomic structure that guide evolutionary studies (^{Platt et al., 2015}). One of the most diverse species groups is the P. mexicanus group (^{Pérez-Consuegra & Vázquez-Domínguez, 2015}, ²⁰¹⁷). Originally, it included only 4 species (^{Osgood, 1909}), but several taxa have since been added (^{Bradley et al., 2007}; ^{Carleton, 1989}; ^{Pérez-Consuegra & Vázquez-Domínguez, 2015}; ^{Rogers & Engstrom, 1992}). In later decades, the P. mexicanus group fluctuated in diversity, ranging from 7 to 17 species (^{Wade, 1999}), but the most recent molecular investigations, based on cyt-b, showed that this group includes at least 15 allopatric clades that are morphometrically discernible (Pérez-Consuegra & Vázquez-Domínguez, 2015, 2017). The P. mexicanus species group currently includes P. mexicanus (from San Luis Potosí to Chiapas, in Mexico), P. bakeri (Guatemala), P. carolpattonae (Chiapas and Guatemala), P. gardneri (Guatemala), P. nicaraguae (Honduras and Nicaragua), P. nudipes (Costa Rica), P. grandis (Guatemala), P. guatemalensis (Guatemala), P. gymnotis (Chiapas and Guatemala), P. salvadorensis (El Salvador, Honduras, and Guatemala), P. tropicalis (Belize and Guatemala), P. zarhynchus (Chiapas), in addition to 3 undescribed clades previously named as “G” (Chiapas), “O” (Guatemala), and “L” (Guatemala and Honduras) (^{Álvarez-Castañeda et al., 2019}; ^{Lorenzo et al., 2016}; ^{Pérez-Consuegra & Vázquez-Domínguez, 2015}, ²⁰¹⁷).

Despite recent advances, discrepancies and information gaps remain. For instance, based on morphological descriptions, ^{Ramírez-Pulido et al. (2005)}, ^{Trujano-Álvarez and Álvarez-Castañeda (2010)}, and ^{Pardiñas et al. (2017)}, recognized 7 subspecies in P. mexicanus (P. m. mexicanus [from San Luis Potosí to Oaxaca], P. m. angelensis [coastal region in southeastern Guerrero and southern Oaxaca], P. m. azulensis [Oaxaca’s mountains at the east of the Isthmus of Tehuantepec], P. m. putlaensis [west-central Oaxaca], P. m. saxatilis [southern Chiapas and Guatemala], P. m. teapensis [southern Veracruz, Tabasco, northern Chiapas], and P. m. totontepecus [northeastern Oaxaca]), but ^{Ramírez-Pulido et al. (2014)} also recognized P. m. tehuantepecus (Oaxaca’s lowlands at the south of the Isthmus of Tehuantepec) (Fig. 1). Additionally, recent works showed that the southernmost distribution of P. mexicanus is northern Chiapas (^{Pérez-Consuegra & Vázquez-Domínguez, 2015}, ²⁰¹⁷), so that P. m. saxatilis could represent an isolated population (its type locality is in Guatemala), but the names P. nicaraguae, P. salvadorensis, and P. tropicalis were resurrected from synonymy with P. m. saxatilis, when these 3 Central American clades were found to be distant relatives of P. mexicanus (^{Pérez-Consuegra & Vázquez-Domínguez, 2015}), meaning that P. m. saxatilis should no longer be considered a subspecies of P. mexicanus. On the other hand, a research based on cranial morphological variation suggested that P. m. teapensis and P. m. totontepecus could be treated as independent species (^{Zaragoza-Quintana, 2005}), and even though mitochondrial DNA studies showed that P. mexicanus is paraphyletic (^{Bradley et al., 2007}; ^{Rogers & Engstrom, 1992}; ^{Wade, 1999}), and that P. m. teapensis is closely related to P. gymnotis and not to P. mexicanus (^{Pérez-Consuegra & Vázquez-Domínguez, 2017}), it was relegated to subspecific status within P. mexicanus (^{Pardiñas et al., 2017}; ^{Pérez-Consuegra & Vázquez-Domínguez, 2015}, ²⁰¹⁷).

Southern Mexico is home to several subspecies of the Mexican endemic deer mouse, Peromyscus mexicanus, that have not been included in previous molecular systematics revisions. Those populations deserve particular attention from systematists because the heterogeneous topography in Southern Mexico has promoted a diverse small mammal fauna in which numerous taxonomic arrangements have been confused (^{Guevara et al., 2014}; ^{Guevara & Cervantes, 2014}; ^{León-Paniagua et al., 2007}; ^{Vallejo & González-Cózatl, 2012}; ^{Woodman & Timm, 1999}). Here, we combined previously published genetic data and new mitochondrial DNA sequences from recently collected P. mexicanus from Oaxaca and Tabasco (P. m. teapensis, P. m. tehuantepecus, P. m. totontepecus) (Fig. 1) to assess systematic relationships and to analyze if mitochondrial DNA supports past morphological, molecular and/or taxonomic hypotheses that could clarify some remaining taxonomic issues in the Mexican deer mouse (^{Bradley et al., 2007}; ^{Pérez-Consuegra & Vázquez-Domínguez, 2017}; ^{Ramírez-Pulido et al., 2014}; ^{Rogers & Engstrom, 1992}; ^{Wade, 1999}; ^{Zaragoza-Quintana, 2005}).

Figure 1 Map of southern Mexico showing the geographic ranges of the Peromyscus mexicanus subspecies, with samples analyzed in this work (yellow circles) and type localities of each subspecies (red diamonds) shown. Map modified from ^{Pardiñas et al. (2017)}, ^{Pérez-Consuegra and Vázquez-Domínguez (2017)}, and ^{Trujano-Álvarez and Álvarez-Castañeda (2010)}. 

Materials and methods

All samples came from specimens hosted in scientific collections, and no ad hoc fieldwork was performed for this research. Voucher specimens are deposited in the following mammal collections: Universidad Nacional Autónoma de México, Facultad de Ciencias, Mexico City, Mexico (MZFC); Louisiana State University, Museum of Natural Science, Baton Rouge, USA (LSUMZ); University of California, Museum of Vertebrate Zoology, Berkeley, USA (MVZ); and Texas Tech University, Museum of Texas Tech University, Lubbock, USA (TTU) (Appendix). We sequenced specimens referred to Peromyscus mexicanus and used morphological traits defined by ^{Álvarez-Castañeda et al. (2015)} to verify their identity. We also followed the original descriptions, and used the geographic origin of samples to identify them at the subspecies level (^{Goodwin, 1956}, ¹⁹⁶⁴; ^{Merriam, 1898}; ^{Osgood, 1904}). We identified 2 specimens from Veracruz as P. m. mexicanus (LSUMZ 36423, MZFC 11166), 2 specimens from Teapa, Tabasco (type locality of P. m. teapensis) as P. m. teapensis (MZFC-PMM068, MZFC- PMM081), 2 specimens from Nizanda, Oaxaca (Oaxaca’s lowlands at the south of the Isthmus of Tehuantepec) as P. m. tehuantepecus (MZFC 9288, MZFC 9308), and 2 specimens from northeastern Oaxaca as P. m. totontepecus (MZFC 8689, MZFC 14085). The missing subspecies not analyzed in this work were P. m. angelensis, P. m. azulensis, and P. m. putlaensis.

We sequenced these 8 specimens, along with 2 P. carolpattonae from Chiapas (MZFC 9736, MZFC 9737), Mexico; 1 P. gardneri from Huehuetenango, Guatemala (MVZ 223293); 1 Peromyscus “L” (sensu ^{Pérez-Consuegra & Vázquez-Domínguez, 2017}) from Santa Rosa, Guatemala (MVZ 235925); 1 Peromyscus “O” (sensu ^{Pérez-Consuegra & Vázquez-Domínguez, 2017}) from El Progreso, Guatemala (MVZ 223386); 2 P. nudipes from Cartago and San José, Costa Rica (LSUMZ 28348, LSUMZ 29483); 1 P. salvadorensis from Chalaltenango, El Salvador (MZFC 10917); and 1 P. tropicalis from Izabal, Guatemala (TTU 62079). Additionally, we sequenced 1 P. megalops from Guerrero, Mexico (MZFC 11449), that was used as an outgroup (^{Bradley et al., 2007}).

We extracted whole genomic DNA using a Qiagen DNEasy Blood & Tissue kit (Qiagen, Germantown, Maryland), following the manufacturer’s recommended protocols. Through polymerase chain reaction (PCR), we amplified the complete mitochondrial gene cytochrome b (cyt-b) using the primers MVZ05 (^{Smith & Patton, 1993}) and H15915 (^{Irwin et al., 1991}). Each PCR had a final reaction volume of 13 μL and contained 6.25 μL of GoTaq Green Master Mix (Promega, Madison, WI, USA), 4.75 μL of H20, 0.5 μL of each primer [10μM], and 1 μL of DNA stock. The PCR thermal profile included 2 minutes of initial denaturation at 95 °C, followed by 38 cycles of 30 seconds of denaturation at 95 °C, 30 seconds of annealing at 50 °C, and 1:08 minutes for the extension at 72 °C. We included a 5-minute final extension step at 72 °C. We visualized 3 μL of each PCR product using electrophoresis in 1% agarose gels, stained with SYBR Safe DNA Gel Stain (Life Technologies, Carlsbad, CA, USA). We cleaned each PCR product using 1 μL of a 20% dilution of ExoSAP-IT (GE Healthcare Bio-Sciences Corp. Piscataway, NJ, USA) incubated for 30 minutes at 37 °C followed by 15 minutes at 80 °C. Samples were cycle-sequenced using 6.1 μL of H20, 1.5 μL of 5x buffer, 1 μL of 10μM primer, 0.4 μL of ABI PRISM Big Dye v. 3.1 (Applied Biosystems, Foster City, CA, USA), and 1 μL of the cleaned template. The cycle-sequencing profile included 1 minute of initial denaturation at 96 °C, followed by 25 cycles of 10 seconds for denaturation at 96 °C, 5 seconds for annealing at 50 °C, and 4 minutes for the extension at 60 °C. Cycle sequencing products were purified using an EtOH-EDTA precipitation protocol and were read with an ABI 3130xl genetic analyzer (Applied Biosystems, Foster City, CA, USA). DNA sequences were edited, aligned, and visually inspected using Mega X (^{Kumar et al., 2018}) and FinchTV 1.4 (^{Patterson et al., 2004}). We added 41 sequences available in GenBank from the P. mexicanus species group (^{Bradley et al., 2007}; ^{Pérez-Consuegra & Vázquez-Domínguez, 2015}, ²⁰¹⁷) (Appendix).

We used maximum likelihood (ML) and Bayesian inference (BI) to estimate phylogenetic relationships of the P. mexicanus species group. We selected the best partition (initially divided by codon position) and the best-fit substitution model in PartitionFinder 2 (^{Lanfear et al., 2016}), through an exhaustive search among all available models in MrBayes 3.2 (^{Ronquist et al., 2012}), and using the Bayesian Information Criterion (BIC). We used this result for parameterizing both ML and BI. In IQ-TREE 1.6.12 (^{Nguyen et al., 2015}), we estimated the ML gene tree, and assessed branch support with 1,000 nonparametric bootstraps. In MrBayes 3.2 we used 3 hot chains and 1 cold chain in 2 independent runs of 10 million generations, sampling every 1,000 iterations. We checked for convergence of Metropolis-Coupled Markov Chain Monte Carlo results by examining trace plots and confirming in Tracer 1.7 (^{Rambaut et al., 2018}) that effective sample sizes for all parameters were > 200. The final topology was obtained by calculating a majority rule consensus tree from the posterior, after removing a burn-in of 25%.

To test alternative topologies, we inferred the best tree consistent with a constraint that forced the monophyly of P. m. teapensis, P. m. tehuantepecus, or P. m. totontepecus, and all other P. m. mexicanus samples (i.e., 3 distinct constraints). These were estimated in MrBayes, as described above. To compare the majority rule consensus trees of unconstrained BI and the best trees consistent with the constraint, we used the Shimodaira-Hasegawa test (^{Shimodaira & Hasegawa, 1999}) as implemented in phangorn (^{Schliep, 2011}). We compared the likelihoods of each topology assuming an HKY+I+G substitution model (calculated in PartitionFinder 2) and 10,000 bootstrap replicates. We performed analyses with and without optimizing the rate matrices, gamma distribution, invariable sites, and base frequencies.

To evaluate levels of differentiation among clades in the P. mexicanus species group, we calculated p-distances in Mega X, using the pairwise deletion option and the Kimura 2-parameter model (^{Kimura, 1980}), and we draw a heat map using the ggplot2 R package (^{Wickham, 2011}). Based on preliminary analyses, we found that some of our samples could represent an independent mitochondrial clade, so we analyzed its levels of genetic structure to understand its distinctness better. Thus, in Arlequin 3.5 (^{Excoffier & Lischer, 2010}), we used an analysis of molecular variance (AMOVA) to estimate genetic diversity among clades (FCT), among localities within clades (FSC), and within localities (FST), taking into account the molecular distance between haplotypes (^{Pérez-Consuegra & Vázquez-Domínguez, 2015}). We used the matrix of pairwise differences corrected with the Kimura 2-parameter model and used 10,000 permutations to evaluate the significance of our results (α = 0.05). Finally, to further characterize genetic diversity, we used DnaSP 5.10 (^{Librado & Rozas, 2009}) to calculate the number of segregating sites, the number of haplotypes, haplotype diversity (Hd), and nucleotide diversity (π) for each clade.

Results

We sequenced 1,143 base pairs of the mitochondrial cyt-b in specimens of Peromyscus megalops (1), P. carolpattonae (2), P. gardneri (1), Peromyscus “L” (sensu^{Pérez-Consuegra & Vázquez-Domínguez, 2017}; 1), Peromyscus “O” (sensu^{Pérez-Consuegra & Vázquez-Domínguez, 2017; 1}), P. m. mexicanus (2), P. m. teapensis (2), P. m. tehuantepecus (2), P. m. totontepecus (2), P. nudipes (2), P. salvadorensis (1), and P. tropicalis (1) (see Appendix for museum catalog numbers and GenBank Accessions). With these and previously published sequences, we analyzed 59 individuals in total. Our alignment was complete at > 92% of positions, contained 354 variable characters and 302 parsimony-informative characters.

The best evolutionary model scheme used was K80+I+G, HKY+I, and HKY+G applied to the first, second, and third codon positions, respectively. Topologies from ML and BI had only slight differences. The ML tree showed an unsupported but resolved clade where P. guatemalensis was sister to P. grandis + P. bakeri, and all of them were sister to P. carolpattonae; however, in the BI topology P. guatemalensis, P. grandis + P. bakeri, and P. carolpattonae formed a polytomy. Similarly, the ML tree showed an unsupported clade where the sample MZFC 14085 (P. m. totontepecus) was sister to P. m. teapensis, and these specimens were sister to P. m. tehuantepecus; however, in the BI the sample MZFC 14085, P. m. teapensis, and P. m. tehuantepecus appeared as another polytomy.

There are 4 main clades within this group: I) P. tropicalis (Guatemala), II) P. nudipes (Costa Rica), III) 11 subclades that mainly inhabit highlands from Chiapas, Mexico to Nicaragua, and IV) 3 subclades that mainly inhabit low and mid altitudes from Mexico to Guatemala (Fig. 2). Even though bootstrap values and posterior probabilities supported the monophyly of these 4 clades, their relationships were highly supported in the BI tree, but unsupported in the ML analysis. However, the sister relationship between clades II and III did not receive significant support values in either ML or BI (Fig. 2). Within clade III, we found that P. zarhynchus was sister to P. gardneri + Peromyscus “G”; P. salvadorensis was sister to Peromyscus “L”; and P. salvadorensis + Peromyscus “L” was sister to an unresolved subclade that includes P. carolpattonae, P. guatemalensis, and P. grandis + P. bakeri. Besides, P. bakeri, P. carolpattonae, P. grandis, P. guatemalensis, Peromyscus “L”, and P. salvadorensis were sister to P. nicaraguae and Peromyscus “O”. Within clade IV, P. mexicanus was paraphyletic. Peromyscus m. mexicanus (Hidalgo and Veracruz) was sister to a well-supported subclade that placed samples of P. m. teapensis from Chiapas and Tabasco (TTU 82759, MZFC-PMM068, MZFC-PMM081), P. m. tehuantepecus from southeastern Oaxaca (MZFC 9288, MZFC 9308), and P. m. totontepecus from northeastern Oaxaca (MZFC 8689, MZFC 14085) as sister to P. gymnotis (Pacific slope in Mexico and Guatemala). The constrained analyses, which forced P. m. teapensis, P. m. tehuantepecus, or P. m. totontepecus to be members of P. mexicanus, generated worse likelihoods in all cases, and the Shimodaira-Hasegawa test rejected each constraint with p-values < 0.001 (Table 1).

Figure 2 Majority rule consensus tree obtained from Bayesian analysis of cytochrome b sequences, in the P. mexicanus species group. Support values are shown as posterior probabilities followed by bootstrap values from a maximum likelihood analysis. Support values < 0.8/80 are not shown. Numbers I to IV correspond to the 4 main clades detected in this species group; nomenclature of clades according to ^{Pérez-Consuegra and Vázquez-Domínguez (2017)}. Note that P. m. teapensis, P. m. tehuantepecus, and P. m. totontepecus (in the orange box), are sister to P. gymnotis, and not to P. m. mexicanus. Label tips show country (CR = Costa Rica, ES = El Salvador, GT = Guatemala, HO = Honduras, MX = Mexico, NI = Nicaragua), states/provinces/departments (cg = Cartago, sj = San José; ct = Chalatenango; av = Alta Verapaz, ch = Chimaltenango, ep = El Progreso, es = Escuintla, gt = Guatemala, hu = Huehuetenango, iz = Izabal, ju = Jutiapa, qe = Quetzaltenango, qc = Quiche, sr = Santa Rosa, su = Suchitepequez; ol = Olancho, le = Lempira; cs = Chiapas, gr = Guerrero, hi = Hidalgo, oa = Oaxaca, ta = Tabasco, ve = Veracruz; et = Estelí), and catalog number of each analyzed specimen. 

In the AMOVAs, the percentage of variation within localities was minimal in both comparisons, so samples of P. m. mexicanus, P. gymnotis, P. m. teapensis, P. m. tehuantepecus, and P. m. totontepecus did not show signs of mitochondrial panmixia. When P. gymnotis, P. m. teapensis, P. m. tehuantepecus, and P. m. totontepecus were treated as 1 group, the most variation (53.02%) was attributed to among-locality diversity. However, when P. gymnotis was considered an independent group from P. m. teapensis + P. m. tehuantepecus + P. m. totontepecus, the most significant variation was detected among clades (63.25%), reflecting substantial differentiation between them. These results support that P. m. mexicanus, P. gymnotis, and P. m. teapensis + P. m. tehuantepecus + P. m. totontepecus are 3 significantly different mitochondrial groups (Table 2). In the following analyses, we considered samples of P. m. teapensis + P. m. tehuantepecus + P. m. totontepecus an independent group.

Peromyscus zarhynchus, P. tropicalis, Peromyscus “G”, P. nudipes, and P. nicaraguae showed the largest interspecific genetic distances (10.82 - 11.47%) (Fig. 3). In contrast, the smallest interspecific distances were detected among the highland taxa of clade III, which ranged between 2.9 and 4.5%. In clade IV, the genetic distance between P. m. mexicanus and P. gymnotis was 8.31 (range = 7.51 - 9.41), and between P. m. mexicanus and P. m. teapensis + P. m. tehuantepecus + P. m. totontepecus it was 7.53 (range = 6.5 - 8.72), whereas between P. gymnotis and P. m. teapensis + P. m. tehuantepecus + P. m. totontepecus the distance was 6.13 (range = 5.15 - 7.35) (Fig. 3). Finally, in P. mexicanus, P. gymnotis, and P. m. teapensis + P. m. tehuantepecus + P. m. totontepecus, we found high haplotype diversity (Hd > 0.95), but within each species all haplotypes were similar (π < 0.024, segregating sites ≤ 56) (Table 3).

Figure 3 Heat map showing genetic distances in the P. mexicanus species group. Interspecific and intraspecific comparison are shown above and below the black line, respectively. Sample-specific labels are shown on both axes, showing the country (CR = Costa Rica, ES = El Salvador, GT = Guatemala, HO = Honduras, MX = Mexico, NI = Nicaragua), states/provinces/departments (cg = Cartago, sj = San José; ct = Chalatenango; av = Alta Verapaz, ch = Chimaltenango, ep = El Progreso, es = Escuintla, gt = Guatemala, hu = Huehuetenango, iz = Izabal, ju = Jutiapa, qe = Quetzaltenango, qc = Quiche, sr = Santa Rosa, su = Suchitepequez; ol = Olancho, le = Lempira; cs = Chiapas, gr = Guerrero, hi = Hidalgo, oa = Oaxaca, ta = Tabasco, ve = Veracruz; et = Estelí), and catalog number of each analyzed specimen. Nomenclature of clades according to ^{Pérez-Consuegra and Vázquez-Domínguez (2017)}. 

Discussion

Relationships between members of the P. mexicanus species group have been confused because of relatively conserved morphology and karyotypes (^{Rogers & Engstrom, 1992}). However, more recent molecular analyses of cyt-b identified 15 clades that subsequently were re-defined as discrete and allopatric populations fragmented by topography (^{Pérez-Consuegra & Vázquez-Domínguez, 2015}, ²⁰¹⁷). Here, we analyzed samples from P. m. teapensis, P. m. tehuantepecus, and P. m. totontepecus, subspecies that have not had their systematic placement tested with molecular data. We detected a previously unknown clade and showed that P. mexicanus is a paraphyletic assemblage. Peromyscus m. teapensis + P. m. tehuantepecus + P. m. totontepecus were sister to P. gymnotis, not to any other P. mexicanus (Fig. 2, Table 1). These 3 clades are highly divergent (Table 2, and genetic distances from 5.15 to 9.41) compared with many other interspecific mitochondrial distances within the P. mexicanus species group (Fig. 3) (^{Álvarez-Castañeda et al., 2019}; ^{Pérez-Consuegra & Vázquez-Domínguez, 2015}). In the P. mexicanus group, smaller interspecific genetic distances (from 3.59 %) are supported by overall size differences, qualitative cranial traits, univariate and multivariate analyses of craniodental morphometric variables, and allopatric distributions (^{Álvarez-Castañeda et al., 2019}; ^{Lorenzo et al., 2016}; ^{Pérez-Consuegra & Vázquez-Domínguez, 2017}). There are 2 possibilities to solve the paraphyly noticed in P. mexicanus. The first one is to assign the populations of P. m. teapensis, P. m. tehuantepecus, and P. m. totontepecus to P. gymnotis. The second one is to recognize the clade P. m. teapensis + P. m. tehuantepecus + P. m. totontepecus as a different species, and we agree with this hypothesis considering their monophyly and high levels of mitochondrial differentiation. Between P. m. teapensis, P. m. tehuantepecus, and P. m. totontepecus (and even including P. m. azulensis, a subspecies not analyzed here, but that is geographically located between P. m. teapensis, P. m. tehuantepecus, and P. m. totontepecus; Fig. 1), P. totontepecus is the oldest available name, so we recommend referring to this clade as P. totontepecus^{Merriam, 1898} and treating P. m. teapensis and P. m. tehuantepecus as a junior synonyms. The status of P. m. azulensis remains unresolved.

Although future research should test our taxonomic hypothesis with phenotypic and more genetic data, previous studies support our findings. Allozyme variation indicated that samples from northwestern Chiapas (P. totontepecus) shared 2 derived allelic character states with P. gymnotis (^{Rogers & Engstrom, 1992}). Mitochondrial sequences of cyt-b, ND3, ND4L, and ND4 also showed a couple of samples from northern Chiapas sister to P. gymnotis, and all of them were sisters to P. mexicanus from Veracruz (^{Wade, 1999}). In the most recent molecular revision of this species group, based on cyt-b, 3 specimens from northern Chiapas and southern Veracruz were treated as P. cf. gymnotis because those individuals were more closely related to P. gymnotis than to P. mexicanus (^{Pérez-Consuegra & Vázquez-Domínguez, 2017}). Moreover, P. m. totontepecus was initially described as a larger and darker subspecies than P. mexicanus (^{Merriam, 1898}). These pieces of evidence reported previously, including allozymes, DNA sequences from 4 mitochondrial genes, and morphological differences, support recognition of P. totontepecus as an independent lineage.

Additionally, there are intraspecific morphological differences in P. totontepecus. For example, P. m. teapensis was considered similar to P. totontepecus but with some modest differences in color and the skull (^{Osgood, 1904}). These morphological differences were also reported in a more recent analysis of cranial variation, but it was suggested that P. teapensis and P. totontepecus should be treated as independent species (^{Zaragoza-Quintana, 2005}). Based on cranial and external measurements, cranial characters, and fur color traits that we used to verify the specimens analyzed here (^{Merriam, 1898}; ^{Osgood, 1904}), we agree with previous reports, which show some intraspecific variation in P. totontepecus. Furthermore, based on their size, color, and geographic range (^{Merriam, 1898}), we clearly differentiated the specimens MZFC 9288 and MZFC 9308 as “P. m. tehuantepecus”, a subspecies not recognized by ^{Ramírez-Pulido et al. (2005)}, ^{Trujano-Álvarez & Álvarez-Castañeda (2010)}, and ^{Pardiñas et al. (2017)}, but that was recently re-considered as a subspecies of P. mexicanus (^{Ramírez-Pulido et al., 2014}). Future studies should analyze interspecific morphological differences between P. mexicanus, P. totontepecus, and P. gymnotis, but also the intraspecific morphological variation detected in P. totontepecus deserves attention (^{Osgood, 1904}; ^{Zaragoza-Quintana, 2005}).

Based on our results, the geographic distribution of P. gymnotis covers the foothills and adjacent coastal plains in southern Chiapas and northern Guatemala, P. mexicanus ranges from San Luis Potosí to central Veracruz (in the Sierra Madre Oriental), and P. totontepecus likely has a continuous distribution from northeastern Oaxaca to northern Chiapas, including southern Veracruz and southern Tabasco (Fig. 4). It was previously reported that the Isthmus of Tehuantepec (lowlands between Chiapas and Oaxaca) has been a significant barrier that deer mice in the P. mexicanus species group have rarely crossed (^{Ordóñez-Garza et al., 2010}). Our findings agree with this hypothesis because only 2 clades are found to the west (P. mexicanus, and P. totontepecus), but 14 clades are restricted to the mountains east of the isthmus (all species in clades I, II, III, in addition to P. gymnotis). These lowlands are well documented as a geographical barrier for other highland mammals such as shrews (^{Woodman & Timm, 1999}), harvest mice (^{Arellano et al., 2005}), crested tailed mice (^{León-Paniagua et al., 2007}), and other deer mice (^{Sullivan et al., 1997}). However, we found an exception when we analyzed samples of P. totontepecus from both sides of the isthmus. Contrary to previous expectations, mitochondrial sequences were quite similar (Fig. 2, Table 2), showing that the Isthmus of Tehuantepec has acted as a geographic barrier in the P. mexicanus species group, but not for all species (see references in ^{Cortés-Rodríguez et al., 2013}). The species that live in the highest elevations in the P. mexicanus group are placed in the clades II and III, and all of them are restricted to the east (^{Pérez-Consuegra & Vázquez- Domínguez, 2015}, ²⁰¹⁷). However, P. mexicanus and P. gymnotis (clade IV) inhabit low and mid-elevations (^{Pérez-Consuegra & Vázquez-Domínguez, 2017}), as is P. totontepecus (collected between 30 to 2,047 m). At least during interglacial periods like the present, it would be logical for the Isthmus of Tehuantepec to represent a lesser barrier to species adapted to low and middle elevations than those adapted to high elevation ecosystems. This ecological difference between species could explain the dual effect of this isthmus on the evolutionary history of these deer mice.

Our taxonomic hypotheses relied on a single-locus dataset. Still, findings of high levels of genetic divergence, backed by previous allozymic and morphological data, in addition to other mitochondrial loci, provide a better picture of the systematics of the P. mexicanus group. Essential aspects of Peromyscus biology, such as their taxonomic and phylogenetic relationships, are still not clear (^{Bradley et al., 2007}; ^{Miller & Engstrom, 2008}; ^{Platt et al., 2015}; ^{Sullivan et al., 2017}), for example, P. m. angelensis, P. m. azulensis, and P. m. putlaensis have not been included in any molecular phylogenies, and their systematic relationships should be analyzed (Fig. 4). Further research will be crucial to understand the number of Peromyscus species and their phylogenetic relationships. Given their commonness and importance as model organisms (^{Bradley et al., 2007}), this should be a high priority.

Figure 4 Geographic ranges and taxonomic changes proposed in this work. We recognize P. gymnotis, P. mexicanus, and P. totontepecus as independent species, but the status and phylogenetic relationships of P. m. angelensis, P. m. azulensis, and P. m. putlaensis are still unclear.