Two new species of Peromyscus (Cricetidae: Neotominae) from the Transverse Volcanic Belt of Mexico

Abstract

Specimens of the Peromyscus boylii species group distributed in the western and northeastern montane regions of Michoacán, México, historically have been assigned to P. levipes. Previous studies indicated that these specimens possessed mitochondrial DNA haplotypes that were distinct from both P. levipes and P. kilpatricki, a recently named species in the P. boylii species group from northeastern Michoacán and western Morelos. Herein karyotypic, DNA sequence, and morphological data were analyzed from those populations to evaluate their taxonomic affinity. Karyotypic data indicated that individuals from western Michoacán (Dos Aguas and Aguililla) and from a newly discovered population in northeastern Michoacán (Zinapécuaro) were chromosomally similar to P. carletoni (FN = 68) but distinct from other taxa assigned to the P. boylii species group. Analyses of cranial characteristics indicated that, relative to other species in the P. boylii species group, two morphologically distinct groups were present that corresponded to the Dos Aguas/Aguililla and Zinapécuaro populations, respectively. The latter population, although represented by a small sample size (n = 5 specimens), appeared to exhibit some trenchant morphological distinctions compared with other cryptic species in the P. boylii group. Phylogenetic analyses (parsimony, Bayesian, and likelihood) of DNA sequences obtained from the mitochondrial cytochrome-b gene indicated that although the individuals from Dos Aguas/Aguililla and Zinapécuaro formed a sister group relationship, they formed monophyletic clades that differed genetically (2.54%)—a level approaching that seen between other sister species of Peromyscus. Further, the Dos Aguas/Aguililla and Zinapécuaro clade was more closely aligned with a clade containing representatives of P. carletoni and P. levipes instead of with those from closer geographic proximities (P. kilpatricki) located in eastern Michoacán. Together, these results indicated that these two populations seemingly represent two undescribed species in the P. boylii species group for which we propose the names Peromyscus greenbaumi for populations in western Michoacán (circa Dos Aguas and Aguililla) and Peromyscus ensinki for populations in northeastern Michoacán (circa Zinapécuaro).

Los especímenes del grupo de especies de Peromyscus boylii distribuidos en las regiones montañosas occidentales y el noreste de Michoacán, México, históricamente fueron asignados a P. levipes. Sin embargo, estudios previos han indicado que estos especímenes poseen haplotipos de ADN mitocondrial que son distintos de P. levipes y P. kilpatricki, una especie recientemente nombrada en el grupo de especies P. boylii del noreste de Michoacán y el oeste de Morelos. Los datos cariotípicos indicaron que los individuos del oeste de Michoacán (Dos Aguas y Aguililla) y de una población recién descubierta en el noreste de Michoacán (Zinapécuaro) eran cromosómicamente similares a P. carletoni (FN = 68), pero distintos de otros taxones asignados al grupo de especies P. boylii. Los análisis de las características craneales indicaron que en relación con otras especies del grupo de especies P. boylii, dos grupos morfológicamente distintos estaban presentes, y que correspondían a las poblaciones de Dos Aguas/Aguililla y Zinapécuaro, respectivamente. Aunque representada por un tamaño de muestra pequeño (n = 5 especímenes), esta última población pareció exhibir algunas diferencias morfológicas en comparación con las otras especies crípticas en el grupo P. boylii. Los análisis filogenéticos (parsimonia, inferencia bayesiana y verosimilitud) de secuencias de ADN obtenidas del gen mitocondrial citocromo-b, indicaron que, aunque los individuos de Dos Aguas/Aguililla y Zinapécuaro tienen una relación de grupo hermano, estos forman clados monofiléticos que difieren genéticamente (2.54%), nivel que se acerca al observado entre otras especies hermanas de Peromyscus. Además, el clado de Dos Aguas/Aguililla y Zinapécuaro está más estrechamente alineado con un clado que contenía representantes de P. carletoni y P. levipes en lugar de aquellos de proximidades geográficas más cercanas (P. kilpatricki) ubicados en el este de Michoacán. Estos resultados indicaron que estas dos poblaciones aparentemente representan dos especies no descritas en el grupo de especies P. boylii para las cuales proponemos los nombres Peromyscus greenbaumi para poblaciones en el oeste de Michoacán (hacia Dos Aguas y Aguililla) y Peromyscus ensinki para poblaciones en el noreste de Michoacán (circa Zinapécuaro).

cryptic species, cytochrome-b gene, karyotype, morphometrics, Peromyscus, P. boylii species group

Especies crípticas, citocromo-b, cariotipo, morfometría, Peromyscus, grupo de especies P. boylii

The Peromyscus boylii species group is a complex of nine named species of deermice found across the southern United States, Mexico, and northern Central America (see Bradley et al. 2017 for most recent synopsis). The P. boylii species group is composed of many cryptic species; many of which have been elevated to species level. Karyotypic diversity in the P. boylii species group (summarized by Houseal et al. 1987) in conjunction with geographic distributions has been instrumental in the taxonomic revisions within previously identified morpho-groups (Schmidly et al. 1988). In particular, the number of autosomal chromosomes in the karyotype (FN) and DNA sequence data from the mitochondrial cytochrome-b gene (Cytb) have become a catalyst for the elevation of multiple taxa to species rank (see Schmidly et al. 1988; Bradley et al. 2000, 2004, 2014, 2017). Bradley et al. (2017) used sequence data from the Cytb gene in conjunction with chromosomal data from karyotypic Group V (Houseal et al. 1987) to describe a new species (Peromyscus kilpatricki) from eastern Michoacán and western Morelos. The recognition of P. kilpatricki assisted in resolving the taxonomy and systematics of the P. boylii species group, especially as it related to populations of P. boylii-like forms from southwestern México. In addition to the description of P. kilpatricki, Bradley et al. (2017) reported a monophyletic clade of DNA sequences from individuals collected in eastern and southwestern Michoacán, México that perhaps represented another undescribed taxon affiliated with the P. boylii species group. Historically, specimens from this region were assigned, under various name combinations, to P. boylii (Osgood 1909; Hooper 1955, 1961, 1968; Schmidly 1973; Schmidly and Schroeter 1974) and later to P. levipes (Carleton 1977, 1979, 1989). The critical work by Carleton (1977, 1979, 1989) and Carleton et al. (1982) transformed the concept of the P. boylii species group and laid the groundwork for deciphering of chromosomal data (see Houseal et al. 1987 and Smith 1990) and ultimately the taxonomic syntheses of Schmidly et al. (1988), Castro-Campillo et al. (1999), and the recent descriptions of new species resulting from the recent DNA sequence data (see below).

This apparently undescribed taxon reported by Bradley et al. (2017), from eastern and southwestern Michoacán, México, possessed a karyotype (karyotypic Group IV, FN = 65–66, 68; Houseal et al. 1987) that was distinct from most other species in the P. boylii species group (beatae, boylii, levipes, kilpatricki, madrensis, schmidlyi, simulus, and stephani). Although, the new karyotype was similar to that reported for individuals of P. carletoni (Bradley et al. 2014) from southern Nayarit, the distinct geographic distribution precluded an association between the two taxa (P. carletoni and the presumably undescribed species) and provided the emphasis for this study. Further, allozymic data (Rennert and Kilpatrick 1987) revealed that the Michoacán populations (corresponding to the FN = 56–P. kilpatricki and FN = 65–66, 68 undescribed taxon) were genetically divergent from each other and from other members of the P. boylii species group. However, no fixed genic differences were observed between these taxa and other species in the P. boylii species group; consequently, the allozyme data were not useful in resolving the taxonomic affinities of these populations relative to each other. Given the diverse morphotectonic provinces in Michoacán, such as the Faja Volcanica Transmexicana and Sierra Madre del Sur (Ferrusquía-Villfranca 1993; Morrone 2019), it is quite possible that these genetically distinct populations are cryptic species within the P. boylii species group.

In this study, efforts were focused on determining the taxonomic status of the P. boylii species group distributed in the high elevation (>1,900 m) pine–oak forests (Pinus spp. and Quercus spp.) of Michoacán and Jalisco, México. Although, Houseal et al. (1987) treated populations from Jalisco, Michoacán, and Nayarit as a single karyotypic group (IV), DNA sequence data (Tiemann-Boege et al. 2000; Bradley et al. 2004) indicated that populations from Ocota and Santa María, Nayarit, and populations from Dos Aguas, Michoacán, were genetically divergent and that the Nayarit samples would prove to be a separate species (P. carletoni, Bradley et al. 2014). Given the similar geographic proximity to other species in the P. boylii species group (beatae, carletoni, kilpatricki, levipes, and schmidlyi) representative samples from throughout this region were included as were previously unreported samples from Michoacán (Aguililla and Zinapécuaro).

Materials and Methods

Samples

Specimens from five naturally occurring populations, four in Michoacán, México (Localities 19, 26, 27, and B, Fig. 1;  Appendix I) and one in Jalisco, México (Locality A, Fig. 1;  Appendix I) form the basis for the DNA sequence and karyotypic components of this study. The first population (Locality 26; Dos Aguas, Michoacán) provided data for the karyotypic (n = 4) and sequence (n = 3) data sets; the second population (Locality 27; Aguililla, Michoacán) provided data for the sequence data set (n = 2); the third population (Locality 19; Zinapécuaro, Michoacán) provided data for the karyotypic (n = 3) and sequence (n = 4) data sets; the fourth (Locality B; Los Reyes, Michoacán; n = 1); and fifth (Locality A; Volcán de Colima, Jalisco; n = 1) populations provided data for the karyotypic data set. In addition, a number of museum specimens, combined with some of those from the sequence and karyotypic data sets, were used in a morphometric analysis of variation among samples in the Mexican state of Michoacán where multiple species of the P. boylii group have been documented (see morphometric data below).

Fig. 1.

Distribution of selected populations and species of the Peromyscus boylii species group occurring in México. Emphasis was placed on depicting the newly described species and its closest phylogenetic allies, taxa with similar karyotypes, and those of relatively close geographic proximity. Closed circles represent collecting localities, numbers refer to samples listed in Appendix I, and letters refer to populations with relevant chromosomal data (obtained from the literature). Localities A, B, C, D, and 26–27 represent the one of the new species (P. species novum1); whereas Locality 19 represents the second new species (P. species novum2).

Given that specimens obtained from Zinapécuaro were collected in 2006 and 2008 and predated Sikes et al. (2016), we followed the methods outlined in the ASM Guidelines (ad hoc Committee on Acceptable Field Methods in Mammalogy 1998; Gannon et al. 2007) and approved by the Texas Tech University Animal Care and Use Committee; all other data were either obtained from the literature in the cases of chromosomal data or were borrowed from natural history museums (tissues for sequencing). Specimen numbers and collection localities are listed in  Appendix I.

Karyotypic data

Three individuals of Peromyscus from Zinapécuaro were karyotyped following methods outlined in Baker and Qumsiyeh (1988). Five to 10 metaphase spreads were examined and photographed for each individual. Karyograms were assembled based on chromosomal morphology presented in Committee for Standardization of Chromosomes of Peromyscus (1977) and Greenbaum et al. (1994) and compared to karyograms and FNs (Table 1; Fig. 2) previously reported by Lee et al. (1972), Schmidly and Schroeter (1974), Houseal et al. (1987), Smith (1990), and Bradley et al. (2014, 2017).

Table 1.

Comparison of karyotypes for members of the Peromyscus boylii species group examined in this study. All chromosomal assessments are based on nondifferentially stained karyotypes as interpreted from comparisons to data presented in or cited by Houseal et al. (1987) and Smith (1990). Only chromosomes that have been identified as biarmed for the P. boylii species group are included. All karyotypes possessed a biarmed condition for chromosomes 1, 22, and 23 (except in some populations of P. beatae—see Davis et al. 1986). Abbreviations are as follows: a = acrocentric, b = biarmed, and p = polymorphic. References: 1 = Lawlor (1971), 2 = Lee et al. (1972), 3 = Schmidly and Schroeter (1974), 4 = Carleton et al. (1982), 5 = Houseal et al. (1987), 6 = Smith (1990), 7 = Bradley et al. (2004), 8 = Bradley et al. (2014), 9 = Bradley et al. (2017), and 10 = this study. *Bradley et al. (2014) discusses an aberrant karyotype for P. carletoni. Three populations of P. kilpatricki are featured due to their geographic proximity to the new species.

Taxon	FN	Chromosome										Reference
		2	3	4	5	6	7	8	9	10	13
P. beatae	48–54	p	a	a	a	a	a	a	a	a	a	5, 6
P. boylii	52	a	a	a	a	a	a	a	a	a	a	5
P. carletoni	66*	b	b	a	b	b	b	a	b	b	a	4, 8
P. kilpatricki (Zitácuaro, Michoacán; n = 1)	56	b	a	a	a	a	a	a	b	a	a	8, 9
P. kilpatricki (Los Azufres, Michoacán; n = 3)	56	b	a	a	a	a	a	a	b	a	a	5
P. kilpatricki (Pátzcuaro, Michoacán; n = 1)	56	b	a	a	a	a	a	a	b	a	a	5
P. levipes	56–60	b	p	a	a	p	a	a	b	a	a	5, 6
P. madrensis	54	a	a	a	a	a	a	a	b	a	a	4
P. schmidlyi	54–56	p	a	a	a	a	a	a	p	a	a	2, 3, 7
P. simulus	52	a	a	a	a	a	a	a	a	a	a	4
P. stephani	52	a	a	a	a	a	a	a	a	a	a	1
P. new species (Dos Aguas, Michoacán; n = 5)	65, 66, 68	b	b	a	b	b	p	a	b	p	p	5
P. new species (Zinapécuaro, Michoacán; n = 4)	68	b	b	a	b	b	b	a	b	b	b	8, 9, 10
P. new species (Volcán de Colima, Jalisco; n = 4)	68	b	b	a	b	b	b	a	b	b	b	5
P. new species (Los Reyes, Michoacán; n = 1)	66	b	b	a	b	b	b	a	b	b	a	3, 5

Fig. 2.

Phylogenetic tree generated using Bayesian methods (MrBayes; Ronquist et al. 2012) and the GTR+I+G model of evolution. Nodal support values obtained from the Bayesian, maximum likelihood, and parsimony analyses are depicted above branches with bootstrap values (>65) obtained from the parsimony analysis placed to the left of the first slash; bootstrap values (>65) for the maximum likelihood analysis placed between the two slashes; and clade probability values (those ≥0.95 are indicated by an asterisk *) and placed to the right of the second slash.

Sequence data

Mitochondrial DNA was isolated from approximately 0.1 g of frozen liver tissue using the DNeasy Blood and Tissue Kit (Qiagen, Valencia, California). The entire Cytb gene (1,143 bp) was amplified using the polymerase chain reaction method (PCR; Saiki et al. 1988) with the following primers: MVZ05 (Smith and Patton 1993) and PERO3′ (Tiemann-Boege et al. 2000) or primer pairs L14724 with CBH3 (Irwin et al. 1991; Palumbi 1996) and F1 with H15915 (Irwin et al. 1991; Whiting et al. 2003). Thermal profiles for PCR were as follows: initial denaturation at 95°C for 2 min, followed by 35 cycles of denaturation at 95°C for 1 min, annealing at 51°C for 1 min, and extension at 72°C for 2 min, with a final extension at 72°C for 7 min. PCR products were purified with an ExoSAP-IT kit (Affymetrix, Santa Clara, California). Primers used to cycle sequence the products included: WDRAT1100, 400R, 700H, and NEO700L (Peppers and Bradley 2000) and 400F (Tiemann-Boege et al. 2000), H15915, L14724, F1, and CB3H. Cycle sequencing was conducted with BigDye Terminator v3.1 (Applied Biosystems, Foster City, California) using the PCR primers as above. Cycle sequencing reactions were subsequently purified using Sephadex (GE Healthcare, Chicago, Illinois) filtration and centrifugation methods, followed by dehydration and suspension in formamide. Purified sequencing products were analyzed on an ABI 3730xl automated sequencer (Biotechnology Resource Center, Institute of Biotechnology, Cornell University, Ithaca, New York). Sequences were aligned using Sequencher 4.1.2 software (Gene Codes, Ann Arbor, Michigan) and chromatograms were proofed to verify all base changes. Cytb sequences generated herein were deposited in GenBank and are listed in  Appendix I.

Phylogenetic and genetic divergence analyses

A neighbor-joining analysis (PAUP*, Swofford 2003) was conducted on a preliminary data set containing 91 individuals from the P. boylii species group: P. beatae, n = 13; P. boylii, n = 17; P. carletoni, n = 14; P. kilpatricki, n = 7; P. levipes, n = 6; P. madrensis, n = 1; P. schmidlyi, n = 22, P. simulus, n = 1; P. sp., n = 9; and P. stephani, n = 1 to verify taxonomic assignment, eliminate duplicate haplotypes, and confirm monophyly of clades and taxa. From this effort, the final data set used in subsequent analyses included 34 individuals: P. beatae, n = 3; P. boylii, n = 5; P. carletoni, n = 2; P. kilpatricki, n = 4; P. levipes, n = 4; P. madrensis, n = 1; P. schmidlyi, n = 3, P. sp., n = 10; P. simulus, n = 1; and P. stephani, n = 1. Peromyscus gratus was selected as the outgroup taxon for all sequence analyses based on the putative phylogenetic relationships of the genus Peromyscus presented in Bradley et al. (2007). In addition, a single representative of each species from the P. aztecus species group (aztecus, cordillerae, evides, hylocetes, spicilegus, and winkelmanni) was included based on their chromosomal similarities or geographic proximity (see Bradley et al. 2000, 2004, 2007, 2014, 2017; Houseal et al. 1987; Smith 1990; Tiemann-Boege et al. 2000).

A parsimony analysis (PAUP*, Swofford 2003) was conducted using equally weighted characters and variable nucleotide positions treated as unordered, discrete characters with four possible states; A, C, G, or T. Phylogenetically uninformative characters were excluded and the heuristic search option in PAUP* (Version 4.0a169; Swofford 2003) was used to identify the most-parsimonious trees. A strict consensus tree was generated from the available trees and the bootstrap analysis (Felsenstein 1985) with 1,000 iterations was used to evaluate nodal support.

Eighty-eight maximum likelihood models were evaluated using jModelTest2 (Darriba et al. 2012) in order to determine the model of DNA evolution best fitting the data. The Akaike information criterion identified HKY+I+G as being the most appropriate model (−lnL 5312.42), relative to complexity of model, for this data set. A likelihood analysis was performed using the program RAxML (Version 8.2.12, Stamatakis 2014) and the following parameters: base frequencies (A = 0.3272, C = 0.2986, G = 0.1091, and T = 0.2652), proportion of invariable sites (I = 0.5970), and gamma distribution (G = 1.1400) and the GTR+I +G because the RAxML program uses only the GTR model. Nodal support was evaluated using the bootstrap method (10,000 iterations, Felsenstein 1985), with bootstrap values (BS) ≥ 65 used to indicate moderate to strong nodal support.

The MrBayes program (Ronquist et al. 2012) was used to examine DNA sequences under a maximum likelihood framework and Bayesian inference methods were used to obtain clade probabilities values indicative of nodal support. The analysis was run with the following options: 4 Markov-chains, 10 million generations, and sample frequency = every 1,000th generation; additional parameters were estimated within the analysis. Following a visual inspection of likelihood scores, the first 1,000,000 trees were discarded and a consensus tree (50% majority rule) was constructed from the remaining trees.

Two methods were used to assess the magnitude of genetic differentiation and to approximate the timing of divergence events. In the first method, the Kimura 2-parameter model of evolution (Kimura 1980) was used to estimate genetic distances among selected clades and taxa, with identical haplotypes removed from the analysis and only one individual per population was included. These values were used to assess levels of genetic divergence relative to the genetic species concept following criteria outlined in Bradley and Baker (2001) and Baker and Bradley (2006). For example, several studies (Bradley et al. 2004, 2014, 2016, 2017) have shown that species of deermice generally differ by genetic divergence values ranging from 2% to 5%. These values were used as a benchmark to evaluate whether levels of genetic divergence between the FN = 65—66, 68 forms and other members of the P. boylii were indicative of an undescribed species.

In the second method, a molecular clock test, performed in MEGA X (Kumar et al. 2018), failed to reject a strict molecular clock; therefore, divergence dates were estimated using BEAST v2.6.3 (Bouckaert et al. 2019). For the BEAST analyses, branch lengths were scaled but the topology was constrained to the ML tree. A prior lognormal distribution was placed on root height and used to constrain the divergence date estimates of the overall tree to ~3.5 MY with a σ value of 0.5. Calibrations were based on the fossil date proposed for P. gratus (Platt et al. 2015) and followed methods outlined in previous studies (Ordóñez-Garza et al. 2014; Platt et al. 2015; Sullivan et al. 2017; Bouckaert et al. 2019). Three test runs including the GTR+I+G model of evolution, 10,000,000 generations, and a 25% burn-in were used to optimize the analysis and to determine final parameters. A final run used the GTR+I+G model of evolution, 50,000,000 generations, and a 50% burn-in. Log and tree files were then combined to generate divergence date estimates and to produce a maximum clade credibility tree. The program Tracer (Bouckaert et al. 2019) was used to examine for sufficient mixing, convergence stability, and effective sample size >200 for all parameters. Finally, the program TreeAnnotator (Bouckaert et al. 2019) was used to obtain an estimate of the phylogenetic tree.

Morphometric data

Eighteen cranial measurements (defined in Carleton et al. 1982) were recorded in millimeters (mm) from adults or were obtained from previous studies (Bradley et al. 2014, 2017). Only adults (age classes IV–VI), identified as such based on patterns of molar tooth wear (Schmidly 1973), were included in this study. Measurements are as follows: greatest length of skull (GLS), length of auditory bulla (LAB), postpalatal length (PPL), length of mesopterygoid fossa (LMF), palatal length (PL), length of incisive foramen (LIF), length of molar toothrow (LMT), greatest zygomatic breadth (ZB), mastoidal breadth (MB), greatest breadth across molars (GBM), postdental palatal breadth (PPB), greatest width of mesopterygoid fossa (WMF), depth of braincase (DB), breadth of braincase (BB), least interorbital width (LIW), rostral breadth (RB), nasal length (NL), and rostral length (RL).

For the 18 cranial characters, measurements were included from specimens collected at the following localities in Michoacán that form the basis of this study: vicinity of Dos Aguas (Locality 26; n = 13) and Sierra Barcalosa (Locality 26; n = 2); vicinity of Aguililla (Locality 27; n = 2); Los Reyes (Locality B; n = 2); Los Azufres, Puerta Garnica, and Zitácuaro (Localities D and 23; n = 9); Opopeo, Patzcuaro, and Quiroga (Locality C; n = 4); vicinity of Uruapan (Locality E; n = 43); and vicinity of Zinapécuaro (Locality 19; n = 5). These populations were compared with other populations delimited by previous allozyme, karyotypic, morphologic, and genetic studies (Houseal et al. 1987; Rennert and Kilpatrick 1987; Schmidly et al. 1988; Bradley et al. 2004, 2014, 2017) to represent samples of P. beatae (n = 32), P. boylii (n = 44), P. carletoni (n = 54), P. kilpatricki (n = 55), P. levipes (n = 85), P. schmidlyi (n = 9), and P. simulus (n = 15). For all analyses, sexes were combined following Schmidly et al. (1988).

Statistical analyses of morphologic data

For descriptive and comparative purposes, means, ranges, and standard errors were calculated for each character and species; for all further analyses, characters were log-transformed (natural log). Specimens with missing measurements were excluded from analyses. Statistical tests were evaluated at P < 0.05 and were performed using PAST (Hammer et al. 2001) or R software version 3.2.1 (R Core Team 2014).

The Shapiro–Wilk normality test (Shapiro and Wilk 1965) was performed to test for normality among the data. Among-group (localities) variation was examined using an ANOVA on the 18 morphological characters. Not all the variables were normally distributed; for those that were, the Mann–Whitney pairwise test (Ryan 1959) was performed to determine if statistically significant differences existed between pairwise comparisons of the various taxa (P < 0.05). A Dunn post hoc pairwise test (Dunn 1964) was used to test for significance among those variables that were not normally distributed (P < 0.05).

To best explain the variation of the data in multivariate character space, a principal component analysis (PCA) was performed comparing all the species and samples included in the study. In addition, a discriminant function analysis (DFA) was used to produce a scatter plot of species and samples along the first and second axes producing maximal and second to maximal separation among all groups derived from multigroup discriminant function. A classification analysis, associated with the DFA, was conducted to determine the extent to which individual specimens could be correctly assigned to their respective sample group. A separate DFA and classification analysis, using only specimens of the boylii group from Michoacán, was conducted to determine the degree of morphological divergence where their ranges approach one another in that state.

Results

Karyotypic data

Karyotypes were obtained from four individuals collected at Zinapécuaro, Michoacán in 2006 and 2008. Additional karyotypes from geographic localities of closely related species were obtained from Houseal et al. (1987) and Bradley et al. (2017) and compared with those generated herein (Table 1). The karyotype of specimens from Zinapécuaro, Michoacán (TK150637, TK148852, and TK148854), possessed a diploid number (2n = 48) and fundamental number (FN = 68) similar to those reported for specimens from Dos Aguas (FN = 65–66, 68), Los Reyes (FN = 66), Michoacán, and Volcán de Colima (FN = 68) by Houseal et al. (1987). The FN = 68 karyotype contained three large pairs of biarmed, five medium pairs of biarmed, and three small pairs of biarmed chromosomes; the polymorphisms associated with the FN = 65 and FN = 66 karyotype is described in Houseal et al. (1987). This karyotype also was similar to that reported for populations in Nayarit, México that are now referable to P. carletoni (Houseal et al. 1987; Bradley et al. 2017).

Sequence data

The preliminary analysis conducted on the DNA sequences from nine individuals from three localities in Michoacán (Aguililla, Dos Aguas, and Zinapécuaro), indicated that these samples formed a monophyletic group in all DNA sequence analyses (parsimony, likelihood, and Bayesian) within the P. boylii species group, as defined by Tiemann-Boege et al. (2000) and Bradley et al. (2004, 2014, 2017). Each analysis is discussed in detail below; however, only the topology obtained from the Bayesian inference analysis is shown (Fig. 2). Further, relationships of some taxa included as reference samples were peripheral to the goals of this study; consequently, not all taxa and clades are discussed in detail.

In the Bayesian inference analysis (Fig. 2), a major division was apparent between members of the P. boylii species group, separating those occupying the montane regions of México (Clade A) from those from the arid portions of the southwestern United States, Mexican plateau, and island inhabitants of the Golfo de California (Clade B). Clade A contained eight major nodes, labeled as Clades C–J, with each clade (except Clade G) receiving statistical support (CPV > 0.95). The nine samples from Michoacán formed a strongly supported monophyletic clade (Clade H, CPV = 1.00) that was comprised of two subclades, one containing the five individuals from Aguililla and Dos Aguas (Clade I, CPV = 1.00) and one containing the four individuals from Zinapécuaro (Clade J, CPV = 1.00). Members of Clade H and Clade G (CPV < 0.95), which contained representatives of P. levipes and P. carletoni, formed a sister clade (Clade E, CPV = 0.96); followed by the addition of P. schmidlyi to form Clade D (CPV = 1.00). Samples of P. beatae and P. kilpatricki were sister (Clade C, CPV = 1.00) which then was sister to Clade D (Clade A, CPV = 1.00). Samples of boylii, madrensis, simulus, and stephani were added in a step-wise fashion to form a monophyletic P. boylii species group (Clade B; CPV = 1.00).

The maximum likelihood analysis also produced a topology (not shown) that was essentially identical to the topology obtained from the Bayesian inference analysis. The only difference between the Bayesian and maximum likelihood analyses involved the placement of the samples representing P. schmidlyi. In the maximum likelihood analysis, P. schmidlyi was placed as the basal most taxon within Clade A, instead of being sister to Clade E. Bootstrap support values obtained from the maximum likelihood analysis were superimposed onto the Bayesian topology (Fig. 2).

For the parsimony analysis, 12 equally, most-parsimonious trees (length = 764, homoplasy index = 0.4895, and consistency index = 0.5105) were retrieved. A majority rule consensus tree was generated (not shown) that was similar in topology to the tree obtained in the Bayesian analysis (Fig. 2); consequently, the bootstrap support values from the parsimony analysis were superimposed onto the Bayesian topology.

Molecular dating

A test of the molecular clock model (Kumar et al. 2018) indicated that the rates of genetic change were indicative of a strict molecular clock. The BEAST program estimated a Yule birth rate of 0.96 [based on a 95% highest posterior density (HPD) that ranged from 0.41 to 1.64]. The mean rate of evolution for Cytb gene (depicted as substitutions per site per million years) was 0.0173 (based on a 95% HPD) and ranged from 0.0082 to 0.0281. Divergence estimates (Fig. 3) indicated that the initial split between the outgroup taxa (P. gratus) and members of the P. boylii species group was 4.65 million years ago (mya). Divergence times indicated that the divergence of the P. aztecus species group from the P. boylii species group occurred approximately 4.41 mya. A major division between the montane forms (Clade A) and the arid species (Clade B) occurred 3.04 mya. Within the montane forms (Clades C–J), divergence dates ranged from 2.04 to 0.80 mya; with the nine samples from Michoacán diverging from the common ancestor (Clade E) with the clade containing P. carletoni and P. levipes (Clade G) approximately 1.09 mya and from each other 0.80 mya (Clade H).

Fig. 3.

Phylogenetic time-calibrated ultrametric tree obtained from a BEAST analysis (Version 2.6.3, Bouckaert et al. 2019) of the mitochondrial cytochrome-b gene data set. Scale bars (gray horizontal rectangles) at nodes represent the 95% highest posterior densities. Vertical lines represent the estimated divergence times in million years ago and are used to reflect approximate divergence times for each node.

Genetic divergence values (Table 2) were estimated using the Kimura 2-parameter model of evolution (Kimura 1980). Genetic diversity among individuals representing the new species from western Michoacán (Aguililla and Dos Aguas) was 0.35%; whereas genetic distances among individuals of the new species from eastern Michoacán (Zinapécuaro) were 0.62%. The two new species differed from each other by a genetic distance of 2.54%. The new species from western Michoacán differed genetically from other members in Clade A (P. beatae, P. carletoni, P. kilpatricki, P. levipes, and P. schmidlyi) by 5.88%, 3.22%, 5.99%, 2.93%, and 3.41%, respectively. The new species from eastern Michoacán differed genetically from other members in Clade A (P. beatae, P. carletoni, P. kilpatricki, P. levipes, and P. schmidlyi) by 6.31%, 3.63%, 6.42%, 3.19%, and 3.93%, respectively. Comparatively, genetic divergence values between other closely related species in the P. boylii species group ranged from 2.68% (P. carletoni and P. levipes) to 8.57% (P. boylii and P. kilpatricki).

Table 2.

Average genetic distances estimated using the Kimura 2-parameter model of evolution (Kimura 1980) for selected comparisons of taxa of Peromyscus. Identical haplotypes were removed from the analysis and only one individual per population was included.

Comparison	Average genetic distance
Within species
P. beatae	1.29%
P. boylii	1.22%
P. carletoni	0.79%
P. kilpatricki	2.58%
P. levipes	1.12%
P. schmidlyi	0.77%
P. new species 1 (Dos Aguas/Aguilla)	0.35%
P. new species 2 (Zinapécuaro)	0.62%
Between P. new species 1 and other species in the P. boylii species group
P. new species 1 (Dos Aguas/Aguilla)—P. beatae	5.88%
P. new species 1 (Dos Aguas/Aguilla)—P. boylii	8.21%
P. new species 1 (Dos Aguas/Aguilla)—P. carletoni	3.22%
P. new species 1 (Dos Aguas/Aguilla)—P. kilpatricki	5.99%
P. new species 1 (Dos Aguas/Aguilla)—P. levipes	2.93%
P. new species 1 (Dos Aguas/Aguilla)—P. madrensis	7.72%
P. new species 1 (Dos Aguas/Aguilla)—P. schmidlyi	3.41%
P. new species 1 (Dos Aguas/Aguilla)—P. simulus	7.51%
P. new species 1 (Dos Aguas/Aguilla)—P. stephani	8.55%
P. new species 1 (Dos Aguas/Aguilla)—P. new species 2 (Zinapécuaro)	2.54%
Between P. new species 2 and other species in the P. boylii species group
P. new species 2 (Zinapécuaro)—P. beatae	6.31%
P. new species 2 (Zinapécuaro)—P. boylii	9.00%
P. new species 2 (Zinapécuaro)—P. carletoni	3.63%
P. new species 2 (Zinapécuaro)—P. kilpatricki	6.42%
P. new species 2 (Zinapécuaro)—P. levipes	3.19%
P. new species 2 (Zinapécuaro)—P. madrensis	7.76%
P. new species 2 (Zinapécuaro)—P. schmidlyi	3.93%
P. new species 2 (Zinapécuaro)—P. simulus	8.45%
P. new species 2 (Zinapécuaro)—P. stephani	9.03%
Between selected species in the Peromyscus boylii species group
P. beatae—P. boylii	8.26%
P. beatae—P. carletoni	5.71%
P. beatae—P. kilpatricki	5.76%
P. beatae—P. levipes	5.59%
P. beatae—P. schmidlyi	5.50%
P. boylii—P. carletoni	7.65%
P. boylii—P. kilpatricki	8.57%
P. boylii—P. levipes	8.48%
P. boylii—P. schmidlyi	7.84%
P. carletoni—P. kilpatricki	5.55%
P. carletoni—P. levipes	2.68%
P. carletoni—P. schmidlyi	3.51%
P. levipes—P. kilpatricki	5.54%
P. levipes—P. schmidlyi	3.08%
P. kilpatricki—P. schmidlyi	6.02%

Morphometric data

Means, ranges, and SEs for the 18 cranial measurements for all species in the P. boylii group are presented in Table 3. ANOVA for the 18 measurements revealed statistically significant differences (P < 0.05) among the taxa in each of the measurements. The Shapiro–Wilk normality test revealed that all but five of the cranial measurements (LN, PPL, ZB, BB, and MB) were normally distributed. Application of the Mann–Whitney test (for normally distributed variables) or Dunn’s post hoc test (for non-normally distributed variables) was used to identify significant differences (P < 0.05) in pairwise comparisons of taxa and samples.

Table 3.

Summary of univariate statistics for Peromyscus samples.

	Mean	Range	SD	Mean	Range	SD	Mean	Range	SD
	P. beatae (n = 32)			P. boylii (n = 44)			P. carletoni (n = 54)
Greatest length of skull	28.37	27.35–29.75	0.59	27.10	25.80–28.35	0.76	27.15	25.40–28.70	0.70
Length of rostrum	11.85	10.10–12.75	0.44	11.01	9.80–12.20	0.51	10.32	9.00–12.00	0.92
Length of nasal	10.61	9.98–11.48	0.38	9.89	8.86–11.20	0.58	10.44	8.86–11.70	0.75
Postpalatal length	9.46	8.80–10.45	0.38	9.08	8.40–9.70	0.37	9.28	8.40–10.00	0.42
Zygomatic breadth	14.18	13.60–15.30	0.38	13.32	12.30–14.95	0.53	13.66	12.70–14.60	0.34
Breadth of braincase	12.87	12.40–13.60	0.27	12.46	11.75–13.35	0.32	12.57	11.80–13.80	0.33
Mastoid breadth	12.04	11.55–12.75	0.28	11.62	11.15–12.50	0.29	11.85	11.00–12.60	0.41
Least interorbital width	4.42	4.20–4.70	0.13	4.32	4.00–4.60	0.15	4.37	4.00–4.70	0.16
Length of molar row	4.41	3.95–4.67	0.14	4.19	3.73–5.38	0.30	4.15	3.73–4.60	0.21
Length of incisive foramen	5.55	5.15–6.10	0.22	4.94	4.15–5.65	0.39	5.32	4.70–5.90	0.30
Length of auditory bulla	5.53	5.13–5.88	0.19	5.30	4.95–5.69	0.18	4.95	3.90–5.69	0.42
Depth of braincase	10.06	9.60–10.45	0.21	9.68	9.10–10.25	0.31	9.66	9.30–10.20	0.23
Length of mesopterygoid fossa	4.90	4.30–5.40	0.24	4.62	4.00–5.20	0.27	4.35	3.80–4.90	0.31
Length of bony palate	4.23	3.90–4.85	0.20	4.20	3.75–4.70	0.22	4.25	3.70–4.80	0.25
Rostral Breadth	4.59	4.20–4.90	0.19	4.55	4.20–4.90	0.19	4.58	4.05–5.20	0.27
Breadth across molars	5.57	5.30–6.00	0.17	5.30	4.45–5.70	0.23	5.34	3.80–5.80	0.30
Postdental palatal breadth	4.11	3.70–4.50	0.20	3.97	3.20–4.35	0.22	3.96	3.30–4.80	0.29
Width of mesopterygoid fossa	2.46	2.20–2.70	0.13	2.33	2.10–4.35	0.15	2.35	1.90–2.80	0.19
	P. ensinki (n = 5)			P. greenbaumi (n = 16)			P. kilpatricki (n = 49)
Greatest length of skull	28.84	28.20–29.75	0.75	28.47	27.75–29.50	0.64	27.15	26.05–28.85	0.63
Length of rostrum	9.46	8.45–10.40	0.80	11.62	10.75–12.35	0.50	10.88	9.60–14.40	0.65
Length of nasal	12.10	11.30–13.10	0.89	10.64	9.60–11.85	0.65	9.75	8.77–11.00	0.43
Postpalatal length	9.56	9.05–10.05	0.41	9.69	9.20–10.30	0.35	9.03	8.20–9.90	0.37
Zygomatic breadth	14.68	14.30–15.00	0.29	14.38	13.65–15.35	0.41	13.58	12.80–14.50	0.36
Breadth of braincase	13.49	13.45–13.50	0.03	13.05	12.45–13.70	0.35	12.67	11.80–14.45	0.42
Mastoid breadth	12.61	12.45–12.90	0.20	12.18	11.25–13.00	0.40	11.59	10.95–12.40	0.30
Least interorbital width	4.51	4.45–4.55	0.05	4.39	4.10–4.70	0.17	4.31	3.95–4.70	0.16
Length of molar row	4.45	4.30–4.65	0.16	4.61	4.29–4.85	0.15	4.39	4.00–4.67	0.14
Length of incisive foramen	6.26	6.10–6.35	0.12	5.47	5.15–5.90	0.21	5.21	4.55–5.80	0.18
Length of auditory bulla	6.84	6.80–6.90	0.05	5.54	5.23–5.79	0.19	5.21	4.75–5.60	0.18
Depth of braincase	10.13	9.70–10.40	0.31	10.05	9.60–10.45	0.26	9.69	9.20–11.20	0.23
Length of mesopterygoid fossa	5.26	5.20–5.35	0.08	4.95	4.35–5.50	0.27	4.85	4.40–5.50	0.26
Length of bony palate	4.54	4.30–4.80	0.21	4.42	4.00–4.90	0.24	4.40	3.65–4.95	0.26
Rostral Breadth	5.03	4.90–5.10	0.09	4.89	4.55–5.75	0.33	4.60	4.15–5.05	0.20
Breadth across molars	5.64	5.55–5.80	0.11	5.63	5.40–6.00	0.19	5.40	5.10–5.70	0.15
Postdental palatal breadth	4.08	4.05–4.10	0.03	3.99	3.75–4.30	0.17	3.96	3.60–4.50	0.17
Width of mesopterygoid fossa	3.05	2.85–3.30	0.19	2.33	2.15–2.50	0.11	2.37	2.10–2.85	0.14
	P. levipes (n = 85)			P. schmidlyi (n = 9)			P. simulus (n = 15)
Greatest length of skull	28.20	25.20–31.15	0.96	26.80	26.30–27.80	0.49	27.12	25.30–28.40	0.78
Length of rostrum	11.37	10.14–13.00	0.57	10.58	9.70–11.30	0.48	10.73	10.00–11.40	0.41
Length of nasal	10.34	8.86–12.65	0.63	9.01	8.40–9.60	0.43	9.47	8.40–10.30	0.50
Postpalatal length	9.58	8.60–11.10	0.48	8.98	8.60–9.40	0.27	9.08	8.00–9.60	0.39
Zygomatic breadth	14.20	13.05–15.60	0.48	13.28	12.60–13.80	0.38	13.89	13.20–14.50	0.41
Breadth of braincase	12.97	12.20–13.85	0.39	12.61	12.40–12.80	0.13	12.35	11.10–11.80	0.24
Mastoid breadth	11.95	11.10–13.15	0.38	12.10	11.70–12.40	0.24	11.46	3.90–4.40	0.20
Least interorbital width	4.45	4.05–4.95	0.18	4.53	4.40–4.80	0.14	4.18	3.60–4.20	0.14
Length of molar row	4.51	3.92–5.13	0.24	4.14	4.00–4.40	0.14	3.83	4.50–5.50	0.16
Length of incisive foramen	5.34	4.35–6.05	0.33	5.54	5.20–5.90	0.25	5.03	5.20–5.90	0.30
Length of auditory bulla	5.39	4.67–5.97	0.24	4.57	4.10–4.90	0.23	5.17	4.10–4.90	0.13
Depth of braincase	9.79	5.40–10.60	0.59	9.37	8.70–9.90	0.39	9.47	8.70–9.90	0.22
Length of mesopterygoid fossa	4.95	4.10–5.90	0.33	3.82	3.50–4.20	0.20	4.80	3.50–4.20	0.24
Length of bony palate	4.44	3.90–5.25	0.26	4.39	4.00–4.70	0.19	4.02	3.70–4.40	0.21
Rostral Breadth	4.75	4.30–5.40	0.23	4.47	4.10–4.80	0.25	4.77	4.50–5.00	0.15
Breadth across molars	5.54	5.05–6.40	0.24	5.26	5.10–5.50	0.13	5.31	4.70–5.90	0.28
Postdental palatal breadth	4.07	3.55–4.50	0.23	3.93	3.80–4.20	0.15	3.95	3.70–4.30	0.15
Width of mesopterygoid fossa	2.36	2.05–2.70	0.12	2.48	2.30–2.70	0.16	2.32	2.10–2.50	0.10

Perusal of the data reveals that with very few exceptions, the external and cranial measurements of the two genetically unique populations from Michoacán (Dos Aguas–Aguililla and Zinapécuaro) average larger than all other species in the P. boylii group, and for many of the cranial measurements the differences are statistically significant (using the pairwise Mann–Whitney or Dunn’s post hoc tests; P < 0.05). The two samples, compared to the other recognized species, are more similar in overall size to P. beatae and P. levipes than the other species in the group.

Pairwise comparisons between the Dos Aguas–Aguililla sample and the other taxa revealed significant differences in 15/18 measurements for the comparison with P. boylii (all but LIW, PPB, and WMF); 14/18 measurements for comparisons with P. schmidlyi (all but MB, LIF, LBP, and PPB) and P. simulus (all but WMF, RB, PPB, and WMF); 13/18 for P. carletoni (all but LN, LIW, LIF, PPB, and WMF). The number of measurements with significant differences were substantially less in the comparisons with P. levipes and P. beatae; only 3 of 18 measurements (MB, LAB, and DB) were significantly different for the comparison with the latter and 4 of 18 (LMT, LBP, RB, and PPB) for the former. In most of the measurements, specimens in the Dos Aguas–Aguililla sample averaged larger than those in the other taxa. The Dos Aguas–Aguililla sample differed significantly from all of the other taxa except for P. levipes in LMR; from all but P. beatae in two measurements (LAB and DB); from all but P. levipes and P. beatae in three measurements (GLS, LR, and BAM); and from all but P. levipes and P. simulus in RB.

Pairwise comparisons of the Zinapécuaro sample with the other taxa revealed an even stronger separation. There were significant differences in all 18 measurements in the comparison with P. simulus; in 17 of 18 measurements in the comparison with P. boylii (all but PPB); in 16 of 18 measurements in the comparisons with P. kilpatricki (all but LMT and LBP) and P. schmidlyi (LIW and LBP); and in 15 of 18 measurements in the comparison with P. carletoni (all but PPL, ZB, and PPB). There were fewer significant differences in the comparisons with the other taxa—only 5 with Peromyscus greenbaumi (LR, LIF, LAB, LMF, and WMF); 9 with P. levipes (LR, LN, BB, MB, LIF, LAB, LMF, RB, and WMF); and 10 in the comparison with P. beatae (LR, LN, MB, LIW, LIF, LAB, LMF, LBP, RB, and WMF). The Zinapécuaro sample differed significantly from all of the other taxa in five measurements (LR, LIF, LAB, LMF, and WMF) and from all but P. greenbaumi in RB and LN.

Although the sample size from Zinapécuaro is small (n = 5), there were trenchant morphometric differences between it and the sample from Dos Aguas–Aguililla (n = 16). With the exception of LR, PPL, and LMR, the Zinapécuaro specimens averaged larger than those from Dos Aguas–Aguililla, and the differences were significantly different in five measurements (LR, LIF, LAB, LMF, and WMF). Furthermore, in four of the five measurements (all but LMF) there was no measurement overlap between these two samples. There also were major differences in the rostrum among specimens of the two samples. The rostrum of the Zinapécuaro specimens made up less than 33% of the length of the skull, and in all five specimens it was significantly shorter than the length of the nasals (mean length of rostrum 9.46 mm vs. mean length of nasals 12.10 mm). The rostrum in all of the Dos Aguas–Aguililla specimens comprised more than 40% of the skull length, and it was always longer than the length of the nasal bones (mean length, 11.62 vs. 10.64 mm). The length of the incisive foramina was another feature with no measurement overlap between the two samples. The range of measurements for the Zinapécuaro sample was 6.1–6.4 mm, whereas in the Dos Aguas–Aguililla sample the range was 5.2–5.9 mm. The same pattern was evident in the width of the mesoptergoid fossa, with Zinapécuaro specimens ranging from 2.8 to 3.3 mm and Dos Aguas–Aguililla ranging from 2.2 to 2.5 mm. The auditory bullae of the Zinapécuaro specimens also were noticeably longer than in the Dos Aguas–Aguililla sample (range of 6.8–6.9 mm in the former compared to 5.3–5.8 in the latter sample).

With the exception of the Zinapécuaro sample, the multivariate analyses depicted less separation among the taxa than did the univariate analysis. The first two components in the PCA accounted for almost 52.2% of the total variation (PC 1 35.6% and PC 2 16.6%). The character loadings were all positive for the first component with the 5 highest being LMF, LAB, LIF, LMR, and LN. For component II, all of the loadings were positive except for four that possessed negative values (GSL, LR, LAB, and DB). The scatter plot for these two components (not shown) revealed broad overlap among the samples except for the Zinapécuaro sample. The Dos Aguas–Aguililla sample was broadly overlapped by the other taxa except for P. carletoni and P. schmidlyi, but it was totally separated in multivariate space from the Zinapécuaro sample.

A DFA was performed to determine the degree of separation in morphometric space among the various taxa and samples and to evaluate how many specimens in each of the samples could be correctly classified. The first two axes of the DFA (see Fig. 4) accounted for 53.4% of the variation among the samples. The loadings for DF 1 were evenly split between positive and negative values (9 vs. 9). Cranial measurements with the highest positive loadings were LAB, LN, LMF, and WMF; those with the highest negative loadings were LR, PPL, DB, and GSL. Loadings for DF 2 also were equally split between positive and negative values. The highest positive loadings were for MB, LIF, GSL, and LAB; the highest negative loadings were for LMR, LR, RB, and BB. The scatter plot of the two axes again showed broad overlap among all of the samples with the exception of the Zinapécuaro sample, which was distinctly separated from the other samples along DF 1. The other samples were broadly overlapped in multivariate space. As in the PCA, there was no overlap among the specimens from Dos Aguas–Aguililla and Zinapécuaro.

Fig. 4.

Plot of the first two discriminant function axes extracted from a discriminant function analysis of 9 species of the Peromyscus boylii species group. This analysis was performed on specimens with complete craniodental measurements. Polygons enclose maximal dispersion of individual specimen scores around centroids for each taxon. Filled circles represent specimens referred to P. beatae, open triangles represent specimens referred to P. boylii, open crosses represent specimens referred to P. carletoni, filled stars represent specimens referred to P. ensinki, open circles represent specimens referred to P. greenbaumi, open stars represent specimens referred to P. kilpatricki, filled squares represent specimens referred to P. levipes, filled triangles represent specimens referred to P. schmidlyi, and filled circles represent specimens referred to P. simulus.

The classification analysis, associated with the DFA, correctly classified 70% of the 309 specimens used in the analysis (Table 4). There were three samples (P. schmidlyi, P. similus, and Zinapécuaro) in which all of the individuals were correctly classified. For other species, 78% of the P. kilpatricki, 74% of the P. carletoni, and 69% of the Dos Aguas–Aguililla specimens (11/16) were correctly classified; four of the latter were misclassified as P. levipes and 1 as P. kilpatricki. 66% of the P. boylii were correctly classified, compared to only 48% of the P. levipes specimens, which was the lowest amount of any of the taxa. Two specimens of P. boylii and 16 P. levipes were misclassified as belonging to the Dos Aguas–Aguililla group. No specimens of the recognized species in the P. boylii group were misclassified into the Zinapécuaro group.

Table 4.

Classification matrix developed from the discriminant function analysis for members of the Peromyscus boylii species group. Numbers reflect where individual specimens were assigned based on analyses of morphometric data.

	beatae	boylii	carletoni	greenbaumi	ensinki	kilpatricki	levipes	schmidlyi	simulus	Total
beatae	29	0	2	0	0	0	1	0	0	32
boylii	4	29	0	2	0	6	3	0	0	44
carletoni	3	7	40	0	0	3	1	0	0	54
greenbaumi	0	0	0	11	0	1	4	0	0	16
ensinki	0	0	0	0	5	0	0	0	0	5
kilpatricki	1	4	1	0	0	38	5	0	0	49
levipes	6	3	2	16	0	17	41	0	0	85
schmidlyi	0	0	0	0	0	0	0	9	0	9
simulus	0	0	0	0	0	0	0	0	15	15
Total	43	43	45	29	5	65	55	9	15	309

Two samples of P. kilpatricki (Uruapan and Zitácuaro), along with those from Zinapécuaro and Dos Aguas–Aguililla, all of which are known from the state of Michoacán, were subjected to a separate DFA to examine the degree of morphological divergence within that state. Specimen plots from these samples are depicted in Fig. 5. The distinctness of the Zinapécuaro specimens stand out in this plot as there is no overlap with specimens from Dos Aguas–Aguililla or those of P. kilpatricki. There is broad overlap among the two P. kilpatricki samples and that from Dos Aguas–Aguililla. The highest character loadings for the first discriminant function in this analysis were LN, LIF, LAB, and WMF. The classification analysis of the 68 Michoacán specimens resulted in 86% correct allocations, including 100% of the Zinapécuaro specimens, 93% of P. kilpatricki (Uruapan sample), 85% of P. kilpatricki (Zitácuaro), and 67% of the Dos Aguas–Aguililla specimens; six of the latter group were incorrectly classified, three each from the two samples of P. kilpatricki.

Fig. 5.

Plot of the first two discriminant function axes extracted from a discriminant function analysis of 3 species of the Peromyscus boylii species group from Michoacán. This analysis was performed on specimens with complete craniodental measurements. Polygons enclose maximal dispersion of individual specimen scores around centroids for each taxon. Filled triangles represent specimens referred to P. ensinki, open circles represent specimens referred to P. kilpatricki (Zitácuaro), open squares represent specimens referred to P. kilpatricki (Uruapan), and closed squares represent specimens referred to P. greenbaumi.

With the exception of the five specimens from Zinapécuaro, it appears that despite a few significant differences among the other taxa in a few univariate comparisons, when variation across all morphometric measurements is considered (multivariate analyses), these differences are overridden by the total morphological variation. This morphometric pattern is somewhat surprising because multivariate approaches typically permit discrimination to a finer degree among similar populations, and they reveal more features of geographic variation than usually seen in one-character comparisons, as exemplified in some of the other cryptic species of the P. boylii group (e.g., see Schmidly et al. 1988; Bradley et al. 2014). Importantly, the character loadings of the cranial measurements in the multivariate analyses were generally consistent with those that proved most useful in the univariate analyses, namely those associated with the rostrum and nasal regions of the skull, along with measurements of the size of the braincase, the length of the incisive foramen and the auditory bullae, and the length and width of the mesopterygoid fossa.

Discussion

Together, results from karyotype data, DNA sequence data, a few cranial measurements, and patterns of geographic distribution suggest that populations (historically referred to P. levipes) from western and northern Michoacán and the Colima-Jalisco border and those from northeastern Michoácan represent two undescribed species of Peromyscus, respectively. Further, based on these results, it appears that a second undescribed species exists among populations of Peromyscus from northeastern Michoacán. These two undescribed taxa from Michoacán (represented by Clades I and J) are included in the montane forms of the P. boylii species group (Fig. 3; Clades C–J), which appear to have diverged approximately 2.04–0.80 mya; their common ancestor (Clade E) with the clade containing P. carletoni and P. levipes (Clade G) approximately 1.09 mya. The two Michoacán lineages then diverged from each other circa 0.80 mya. Below, these populations are formally described as new species.

Peromyscus greenbaumi, new species

Holotype

Texas A&M Biodiversity Research and Teaching Collections (TCWC45304); adult male; skin, skull, and skeleton. Original number NSF 83-782; GK3282 identifies tissue samples deposited in the Texas A&M Biodiversity Research and Teaching Collections, Texas A&M University.

Type locality

México: Michoacán; 11.8 km WSW Dos Aguas, 2,408 m; collected 26 July 1983.

Paratypes

Two males (TCWC45305, GK3289 and TCWC45307, GK3305) and two females (TCWC45309, GK3307 and TCWC45310, GK3308) deposited in the Texas A&M Biodiversity Research and Teaching Collections, Texas A&M University.

Diagnosis

A species of Peromyscus characterized externally by medium size for the genus (mean total length of 10 adults = 207.9 mm); tail about as long as the head and body (for 10 and 15 specimens, respectively, mean of tail = 103.5 mm; mean of head and body = 104.4 mm); hind foot medium (mean = 23.4 mm for 15 adults); and ear medium (mean = 21.67 for 15 adults). According to the color nomenclature of Ridgway (1912), dorsal coloration Clove Brown at tips and Blackish or Mouse Gray at base (color nomenclatural following Ridgway 1912); sides Dresden Brown; venter pelage Dusky Neutral Gray at base and White at tips; feet with Fuscous strip extending from ankle to mid-medapodials; toes White; tail slightly bicolored, Fuscous Black above and Light Grayish-Olive below, scantily haired at base and slightly tufted at tip; ears Raw Umber; and vibrissae Black.

Peromyscus greenbaumi is characterized cranially (based on an adult sample of 16 specimens) by an elongate skull that is almost twice as long as wide (mean greatest length = 28.47 mm; zygomatic breadth = 14.38 mm); a rostrum (mean length = 11.62 mm) that averages 41% of skull length; nasals (mean length = 10.64 mm) that are about 92% of rostral length; a large, rounded braincase (breadth of braincase = 13.05 mm); length of molar tooth row (mean length = 4.61 mm) about 16% of skull length; incisive foramina long (mean length = 5.47 mm) and about 20% of skull length; interorbital constriction with a smooth and not an angular outline; zygomatic arches that are nearly parallel; and large auditory bullae for the P. boylii species group (mean length = 5.54). The karyotype of P. greenbaumi is characterized by a 2n = 48 with 8 pairs of biarmed chromosomes and three pairs that are polymorphic, resulting in an FN of 65, 66, or 68 (see Houseal et al. 1987 for illustrations of this polymorphic karyotype). Distinct mtDNA sequences (Cytb gene) were obtained for P. greenbaumi (GenBank Accession numbers: AF155409, MW684863, MW684864, KX523185, and MW684865).

Distribution

Three species of the P. boylii group occur in Michoacán and their distributions are depicted in Fig. 1. Peromyscus greenbaumi occurs in the Sierra de Coalcomán of western Michoacán, México, beginning in southwestern Michoacán and extended northward into Jalisco along the Colima-Jalisco border, where it occupies the pine–oak forests at elevations from 1,600 to 2,500 m. Peromyscus kilpatricki occupies the highland pine–oak forests of the Purépecha Forest in central Michoacán, and its western limits correspond with the Rio Tepalcatepec Basin in western Michoacán. These two species have been taken at the same locality, 8.4 mi N Los Reyes, in northwestern Michoacán in an area of mesquite-grassland, oak and pine forest (see Remarks below). Peromyscus greenbaumi is separated from P. carletoni to the north by the Rio Lerma Basin and from P. levipes and P. ensinki (described below) to the east by the Valley of Mexico.

Measurements

External measurements of the holotype as taken in the field (in mm) are: total length, 200; tail length, 98; hindfoot, 21; and ear, 20. Cranial measurements were obtained using dial calipers (in mm), accurate at 0.5 mm, and are as follows: GLS, 27.8; ZB, 14.0; LAB, 5.2; PPL, 9.3; LMF, 4.9; LBP, 4.3; LIF, 5.4; LMR, 4.5; MB, 12.0; BAM, 5.5; PPB, 4.3; WMF, 2.3; DB, 9.8; BB, 13.1; LIW, 4.7; RB, 4.7; NL, 9.8; and LR, 10.8. Mean measurements, ranges, and standard errors for additional specimens are presented in Table 3.

Comparisons

A species in the P. boylii species group, resembling the other eight mainland species in external and cranial size and coloration; however, distribution (allopatric), karyotype, and DNA sequence divergence (Cytb gene) preclude confusion. The cryptic nature of the various species in the P. boylii group is confirmed by the morphological analysis.

External measurements average larger than all other species in the P. boylii group except for P. ensinki, which averages larger in total length, tail length, and body length, and P. beatae which is slightly larger in all external measurements than P. greenbaumi. Of the other species, only P. levipes approaches P. greenbaumi in external measurements, although it is slightly smaller in all four of the measurements.

In cranial measurements (see Table 3), compared to P. boylii, P. simulus, and P. carletoni, P. greenbaumi is larger in almost all cranial measurements and the differences are statistically significant (using the pairwise Mann–Whitney or Dunn’s post hoc tests; P < 0.05) for all but three measurements (LIW, PPB, and WMF) in the P. boylii comparison, four measurements in the P. simulus comparison (LMF, RB, PPB, and WMF) and five in the comparison with P. carletoni (LN, LIW, LIF, PPB, and WMF). From P. schmidlyi, P. greenbaumi is larger in all cranial measurements except for LIW, LIF, LBP, and WMF, and the differences are statistically significant for all but four measurements (MB, LIF, LBP, and PPB). In its cranial measurements, P. greenbaumi is most similar to P. levipes, P. beatae, and the newly described P. ensinki. From P. levipes, P. greenbaumi averages larger in all cranial characters but LIW, PPB, and WMF, but the differences are only statistically significant for three measurements (MB, LAB, and DB). From P. beatae, P. greenbaumi averages larger in all cranial measurements except LR, LIW, LIF, DB, PDB, and WMF, but there are only four measurements (LMT, LBP, RB, and PPB) in which the differences are significantly different.

In the multivariate analysis, P. greenbaumi is most distinct from P. schmidlyi, P. simulus, and P. ensinki (see Fig. 4). There is very little morphological differentiation in multivariate space between P. greenbaumi and P. boylii, P. kilpatricki, P. beatae, and P. levipes. The classification analysis revealed five misclassifications of P. greenbaumi (n = 16) with one specimen classified as P. kilpatricki and four specimens as P. levipes.

Peromyscus greenbaumi differs genetically, based on Cytb sequences and by the fundamental number (FN) of the karyotype, from other members of the P. boylii species group. The undescribed taxon differs from P. levipes, to which it was considered conspecific, by a genetic divergence value of 2.93% and it differs from P. ensinki, to which it is a sister species, by a genetic divergence value of 2.54%. A karyotype, based on system the described by the Committee for Standardization of Chromosomes of Peromyscus (1977), depicted a 2n of 48 with eight pairs of biarmed and three pairs of polymorphic autosomes (FN = 65, 66, and 68). Chromosomally, P. greenbaumi (FN = 65, 66, and 68) and P. ensinki (FN = 68) share a similar karyotype, although P. greenbaumi is polymorphic for three medium-sized chromosomes (Table 1).

Etymology

This species is named in honor of Dr. Ira Greenbaum (retired professor, Texas A&M University) who was one of the Principal Investigators on the NSF grant “Chromosomal and biochemical differentiation in Peromyscus: systematics and evolution of the P. boylii group (1981–1984)” that produced much of the karyotypic data, tissue samples, and voucher specimens used in subsequent papers to resolve the systematics of this complex peromyscine group. Dr. Greenbaum also made many other contributions to systematic and evolutionary studies of Peromyscus species.

Habitat

Found in mesic pine/oak forest (Quercus spp. and Pinus spp.) habitat at elevations greater than 1,600 m. This species typically was associated with rock outcroppings, fallen logs, and moist soils and was collected sympatrically with Reithrodontomys sumichrasti, Neotoma mexicana, P. megalops, and P. winkelmanni at the type locality.

Remarks

There is no discrete morphological measurement that will distinguish P. greenbaumi from other species in the P. boylii group, with the exception of P. ensinki, from which it differs as described in the account of the latter. Compared to the other P. boylii taxa, P. greenbaumi is significantly larger (P < 0.05) than all but one of them in three cranial measurements (LAB and DB in P. beatae and LMT in P. levipes). Otherwise, there is some overlap in cranial measurements between P. greenbaumi with two or more of the other species. In both the PCA and DFA of the cranial measurements, P. greenbaumi is broadly included within the character space of all the other taxa except P. schmidlyi and P. simulus. 69% of the specimens of P. greenbaumi (11/16) were correctly identified in the classification analysis. The only reliable way to identify with certainty specimens of P. greenbaumi is to rely on chromosomal or molecular genetic characters. These observations clearly confirm the cryptic nature of the various species in the P. boylii group.

There is a locality in northwestern Michoacán (8.4 mi N Los Reyes) where single individuals with the FN = 66 and FN = 54 karyotypes have been documented together, suggesting that P. greenbaumi and P. kilpatricki possibly occur sympatrically in this area without interbreeding. To test the morphological efficacy of this finding, we included these two specimens in separate PCA and DFA studies with other Michoacán samples. A 2-dimensional plot of the character space in the PCA (not shown), revealed these two specimens to be in opposite quadrants of the graph, with the FN = 66 specimen (the P. greenbaumi karyotype; TCWC 27265) positioned within the character space of the Dos Aguas–Aguililla sample of that species, and the specimen with the FN = 54, P. kilpatricki karyotype (TCWC 27266), positioned just outside of a polygon surrounding specimens of that species from Zitácuaro. The position of the two specimens in the DFA essentially was the same with the exception that the FN = 54 specimen was near samples of P. kilpatricki from both Zitácuaro and Uruapan. A major distinction between the two was that the FN = 54 specimen has a much larger auditory bullae (6.35 mm) than does the FN = 66 specimen (5.05 mm), but overall, the two specimens aligned with the samples of their presumed identity based on karyotypic data which confirms that the ranges of the two species abut in this region of Michoacán. Additional collecting efforts are needed to better refine the distribution of this taxon relative to the other cryptic species in the P. boylii group.

When the holotype and other specimens were collected in the 1980’s, Dos Aguas and other regions in western Michoacán were covered by extensive pine–oak forests. Unfortunately, in recent decades, large tracks of those forests have been converted to pastures, avocado plantations, and some other crops; consequently the landscape has changed dramatically. The impacts of those changes in the populations of P. greenbaum are unknown, but if the species habitat requirements are narrow, it may be endangered and faced with extinction.

Peromyscus ensinki, new species

Holotype

Museum of Texas Tech University MoTTU, (TTU110120); adult male; skin, skull, and skeleton. Original number Robert D. Bradley 2553; TK148853 identifies tissue samples deposited in Natural Sciences Research Laboratory, MoTTU.

Type locality

México: Michoacán; 3.5 km S, 4.8 km E, Zinapécuaro (Santa Cruz), UTM 14Q-311971-2194257, 2,012 m, collected 25 July 2008.

Paratypes

One male (TTU104811, TK150637) and two females (TTU110119, TK148852; TTU110121, TK148854) deposited in the Museum of Texas Tech University.

Diagnosis

A species of Peromyscus characterized externally by medium size for the genus (mean total length = 213.0 mm for three adults); tail about as long as the head and body (mean of tail = 108.0 mm for three adults; mean of head and body = 108.6 mm for four adults); hind foot medium (mean = 22.4 mm for four adults); and ear medium (mean = 20.8 for four adults). According to the color nomenclature of Ridgway (1912), dorsal coloration Sepia at tips and Blackish or Mouse Gray at base (color nomenclatural following Ridgway 1912); sides Cinnamon-rufous; venter pelage Dark Neutral Gray at base and White at tips; feet with Dusky Drab strip extending from ankle to mid-medapodials; toes White; tail slightly bicolored, Fuscous above and Pale Smoke Gray below, scantily haired at base and slightly tufted at tip; ears Deep Mouse Gray; and vibrissae Black.

External measurements of the holotype as taken in the field (in mm) are: total length, 222; tail length, 119; hindfoot, 22; and ear, 21. Cranial measurements were obtained using dial calipers (in mm), accurate at 0.5 mm, and are as follows: GSL, 29.2; ZB, 15.0; LAB, 6.8; PPL, 9.7; LMF, 5.3; LBP, 4.5; LIF, 6.4; LMR, 4.7; MB,12.5; BAM, 5.8; PPB, 4.1; WMF, 3.1; DB, 10.4; BB, 13.5; LIW, 4.6; RB, 5.1; NL, 11.4; and LR, 9.5. Mean measurements, ranges, and standard errors for additional specimens are presented in Table 3.

Peromyscus ensinki is characterized cranially (based on an adult sample of five specimens; see Table 3) by an elongate skull that is twice as long as wide (mean greatest length of skull = 29.75 mm; mean zygomatic breadth = 14.68 mm); a rostrum (mean length = 9.46 mm) that is slightly less than one-third (mean 31.8%) of the skull length; nasal bones that are always larger (mean length = 12.10 mm) than the length of the rostrum; length of molar tooth row (mean = 4.45 mm) about 15% of skull length; long incisive foramina (mean length = 6.26 mm) about 21% of skull length;, very long auditory bullae (mean length = 6.84 mm); a long and wide mesopterygoid fossa (mean length = 5.26 mm; mean width = 3.05 mm); a large, rounded braincase; interorbital constriction with a smooth and not an angular outline; zygomatic arches that are nearly parallel; a karyotype, based on system the described by the Committee for Standardization of Chromosomes of Peromyscus (1977), with a 2n = 48 and 11 pairs of biarmed autosomes (FN = 68). Distinct mtDNA sequences (Cytb gene) were obtained for P. ensinki (GenBank Accession numbers: KF201673, MW684866, and MW684867).

Distribution

Peromyscus ensinki currently is known only from the type locality (see above) and a closely adjacent location (4.8 km S, 3.6 km E Zinapécuaro), both of which are near Lake Zitácuaro in northeastern Michoacán. Specimens of P. kilpatricki have been recorded nearby, including two localities (Los Azufres and Puerto Garnica, Michoacán) approximately 15 km south of Zinapécuaro, and another locality (Quiroga, Michoacán) that is a few km west of the type locality of P. ensinki. No major geographic barriers are evident that would separate the two taxa in this part of Michoacán. Peromyscus ensinki is separated from P. carletoni to the north by the Río Lerma Basin and from P. levipes to the east by the Valley of Mexico. Peromyscus ensinki is separated from P. greenbaumi by the central mountain mass of the Purépecha Highland Forest and the extensive lowland valley of the Río Tepalcatepec Basin that separates the highland forests of central Michoacán from the Sierra de Coalcomán in the western part of the state. Additional collecting efforts are needed to better refine the distribution of this taxon relative to the other species of the P. boylii group in Michoacán.

Comparisons

A species in the P. boylii species group, resembling the other eight mainland species in external and cranial size and coloration; however, distribution (allopatric), karyotype, and DNA sequence divergence (Cytb gene) preclude confusion. The cryptic nature of the various species in the P. boylii group is confirmed by the morphological analysis.

In its external measurements, P. ensinki averages larger in total length, tail and body length, compared to all of the other species in the P. boylii group. Of the other species, only P. greenbaumi, P. beatae, and P. levipes approaches P. ensinki in external measurements. Peromyscus ensinki is substantially larger in its external measurements than P. simulus, P. schmidlyi, P. kilpatricki, P. carletoni, and P. boylii (listed in order of size).

In its cranial measurements and features, P. ensinki is the most distinctive of the species in the P. boylii group. From P. levipes, it is larger in all cranial characters but LR and the differences are statistically significant for all cranial measurements except PPB. From P. beatae, it is larger in all cranial measurements but LR and PPB and in 10 of the measurements (LR, LN, MB, LIW, LIF, LAB, LMF, LBP, RB, and WMF) the differences are statistically significant. Compared to P. boylii, P. simulus, and P. carletoni, P. ensinki is larger in all cranial measurements except for LR, and the differences are statistically significant (using the pairwise Mann–Whitney or Dunn tests; P < 0.05) for all measurements in the P. simulus comparison, all but one measurement (PPB) in the P. boylii comparison, and all but three measurements in the P. carletoni comparison (PPL, ZB, and PPB). Compared to P. schmidlyi, P. ensinki is larger in all cranial measurements except for LR and LIW, and the differences are statistically significant for all but two measurements (LIW and LBP).

In the multivariate analysis, P. ensinki is the most distinctive among all taxa in the P. boylii group (see Figs. 4 and 5). There is no morphological overlap between it and P. greenbaumi and P. kilpatricki, the other two species of the P. boylii group that occur in Michoacán. The classification analysis correctly assigned all five specimens of P. ensinki, and no specimens from the other eight species in the P. boylii group were misclassified as P. ensinki.

Peromyscus ensinki differs genetically, based on Cytb sequences and by the fundamental number (FN) of the karyotype, from other members of the P. boylii species group. The undescribed taxon differs from P. levipes, to which it was considered conspecific, by a genetic divergence value of 3.19% and it differs from P. greenbaumi, to which it is sister, by a genetic divergence value of 2.54%. The karyotype is characterized by a 2n = 48 with 11 pairs of biarmed autosomes (FN = 68). Chromosomally, P. ensinki (FN = 68) share a similar karyotype with P. greenbaumi (FN = 65, 66, and 68) although P. greenbaumi is polymorphic for three medium-sized chromosomes (Table 1).

Etymology

This subspecies is named in honor of Mr. Jan Ensink (deceased M.S. student, Texas A&M University) who was one of the primary members of the field crews (1982–1984) that obtained much of the karyotypic data, tissue samples, and voucher specimens used in subsequent papers to resolve the systematics of the P. boylii species group. Jan was fluent in at least five languages and acted as our interpreter on many field trips to Mexico and Honduras. He was an amazing cook and more often than naught concocted an amazing meal on a Coleman stove from an assortment of ingredients scrounged from the woods and local tienda! More than that, Jan was an exceptional friend whose life ended much sooner than it should have.

Habitat

The two localities where P. ensinki has been found occur in mesic pine/oak forest (Quercus spp. and Pinus spp.) habitat at an elevation of 2,012 m. The specimens were obtained around rock outcroppings, fallen logs, and in moist soils. At the type locality, ensinki was collected sympatrically with Reithrodontomys microdon.

Remarks

Although the sample size for P. ensinki is small (only five specimens), it appears to exhibit some unique features in its cranial architecture and measurements. As mentioned in the description above, it has a small rostrum and long nasal bones compared to the skull length. This, in conjunction with its longer auditory bullae, longer incisive foramen, and wider mesopterygoid fossa, provide unique morphological features that distinguish it from the other species in the boylii group from Michoacán. However, it would be helpful if these morphological differences could be confirmed by a larger sample of specimens.

One of the five specimens of P. ensinki (TTU 110118) is somewhat distinct in its external features, exhibiting a slightly fulvous underbelly compared to the other four specimens. However, in its molecular genetic profile and morphological attributes, this specimen is like the other specimens from Zinapécuaro. A larger sample size will be required to assess the extent to which the fulvous venter might be characteristic of the species.

Nomenclatural statement: A Life Science Identifier (LSID) number was obtained for this publication: urn:lsid:zoobank.org:pub: 07773773-CBD0-448D-A2BA-958AB62EE452

Appendix I

Specimens examined in this study. For each specimen, museum catalogue number (abbreviations for museum acronyms follow Hafner et al. 1997) and GenBank accession number are provided in parentheses. For selected specimens, a collection locality is provided and cross-referenced to Fig. 1. Abbreviations are as follows: Monte L. Bean Life Science Museum (BYU), Mammal Collection, CIIDIR Unidad Durango, Instituto Politécnico Nacional, Durango, Mexico (CRD), Museum of Texas Tech University (TTU), Texas Cooperative Wildlife Collection (TCWC), Universidad Nacional Autónoma de México (UNAM), University of Michigan, Museum of Zoology (UMMZ), and Zadock-Thompson Natural History Collection (ZTNHC). If museum catalogue numbers were unavailable, specimens were referenced with a corresponding collector or other special number (e.g., TK number special number of the Museum of Texas Tech University and CWK special number of Zadock-Thompson Natural History Collection). Several GenBank sequences were not generated in this study and were deposited by Bradley et al. (2000, 2004, 2007, 2014, 2017), Sullivan et al. (1997), and Tiemann-Boege et al. (2000). Localities corresponding to Fig. 1 are provided (in parentheses) are provided only for select members of the P. boylii species group. All specimens utilized in the initial genetic analyses are listed below; specimen selected for the final analyses are indicated by an asterisk (*).

Specimens examined in the DNA sequencing portion of this study

Peromyscus aztecus.—MÉXICO: Veracruz; 8.8 km N Huatusco (TCWC47976, U89968*).

Peromyscus beatae.—EL SALVADOR: Chalatenango, 2.1 mi S La Palma 980 m (ZTNHC: JGO4022, AF131915). GUATEMALA: Sololá, Panajachel (ROM99800, AF131918); Huehuetenango; 10 km NW Sta. Eulalia (ROM98290, AF131919). HONDURAS: Francisco Morazán; 3.2 km NE El Hatillo (TCWC52288, AF131914*). MÉXICO: Chiapas; 18 km S Frontera Comalapa (ROM97642, AF131916); 12 km SE Ixtapa (Locality 34; ROM97582, AF131917). Guerrero; 6.4 km SW Filo de Caballo (Locality 31; TCWC45222, AF131922*). Oaxaca, 6.4 km E Juquila (Locality 32; TCWC45324, AF131920); 3 mi S Suchixtepec (Locality 33; TCWC45385, AF131923*). Veracruz; Xometla (Locality 30; TCWC48060, AF131921); 6.7 km NE, 8.6 km SE Perote (Locality 25; TK150106/TTU105037, MW684859).

Peromyscus boylii.—MÉXICO: Aguascalientes, 6 mi W Rincon de Romos (Locality 17; TCWC48438, AF131924*). Chihuahua, 3.7 mi SW Santa Barbara (Locality 4; TCWC45516, AF131925). Durango; 35.5 km W Cd. Durango (Locality 10; ZTNHC: CWK 1965/AF155413); 15 km N Las Herreras (Locality 6; TK48590, AY322505; TTU75574, AY322506). Jalisco; 2 km NW El Mesconcitos (Locality 18; TTU82688, AY322504*); 30 m W Huejuquilla del Alto (Locality 13; TTU81702, AF155388). Sonora; Isla San Pedro Nolasco (Locality 2; UMMZ117347, AF155387*); 3 km E Yecura, Colegio Yecura (Locality 3; TK148269/TTU110281, MW684860). UNITED STATES: Arizona; Granite Dens Ranch (TK72944/TTU88208, MW684861). California; Orange County, Santa Anna Mountains, Holy Jim Canyon (TTU81248, AF155391*); Monterrey County, Hastings Natural History Reservation (MVZ K. Nutt 120, AF155386); Riverside County, San Bernardino National Forest (TK83719/TTU81355, MW684862). Texas: Culberson County, Sierra Diablo Wildlife Management Area (TTU75833, AF155390*); Jeff Davis County, Mt Livermore Preserve (TTU116513, AF155389). Utah; Garfield County, Henry Mountains, Mt. Pennell, Sidehill Spring, 8,820 ft (MSB123149, AF155392); Washington County; Beaver Dam Wash, 37°00′30″N, 114°07′00″W (TK24389, DQ000478).

Peromyscus carletoni.—MÉXICO: Nayarit; Ocota de la Sierra, 21°50′ N, 104°13′ W (Locality 15; TCWC45206, KF201659*; TCWC45207, AF155410); 70 km N Santa María del Oro, UTM 13Q-559922-2395306 (Locality 16; TTU110264, KF201660; TTU110262, KF201661; TTU110122, KF201662; TTU110124, KF201663; TTU110270, KF201664; TTU110261, KF201665; TTU110263, KF201666; TTU110125, KF201667; TTU110260, KF201668*; TTU110267, KF201669; TTU110265, KF201670; TTU110266, KF201671).

Peromyscus cordillerae.—GUATEMALA: Alta Verapaz, Yalijux Mountain, Chelemhá Reserve, 15°23′09″N, 90°03′44″W, 2,090 m (USNM569872, KF201657*).

Peromyscus ensinki.—MÉXICO: Michoacán; 4.8 km S, 3.6 km E Zinapécuaro (Santa Cruz), 14Q-311982-2194265, 2,012 m (Locality 19; TTU104811, MW684868*); 3.5 km S, 4.8 km E Zinapécuaro (Santa Cruz), 14Q-311971-2194257, 2,012 m (Locality 19; TTU110119, KF201673*; TTU110120, MW684866*; TTU110121, MW684867*).

Peromyscus evides.—MÉXICO: Guerrero; 6.4 km SSW Filo de Caballo (TTU82696, FJ214685*).

Peromyscus gratus.—MÉXICO: Michoacán; Las Minas, 3 KM SW Tuxpan (Catalog number not available TK47810, KF201656*).

Peromyscus greenbaumi.—MÉXICO: Michoacán; 11.8 km WSW Dos Aguas (Locality 26; TCWC45304, AF155409*); 5.0 km SE Dos Aguas, 1,950 m, 180°46.453′N, 102°53.641′W (Locality 26; Catalog number not available—TK45877, MW684863*); 1 km NE Dos Aguas, 2,250 m, 18°48.023′ N, 102°56.081′ W (Locality 26; Catalog number not available—TK45887, MW684864*); 3 km NW Aguililla, 1,780 m, 18°46.238′ N, 102°45.747′ W (Locality 27; Catalog number not available—TK45857, KX523185*); 3.29 km NW Aguililla, 1,900 m, 18°45.351′N and 102°48.585′W (Locality 27; Catalog number not available—TK45865, MW684865*).

Peromyscus hylocetes.—MÉXICO: Michoacán; Estación Cerro Burro, Microondas; 3,270 m (UNAM—catalog number unavailable—TK45309, DQ000481*).

Peromyscus kilpatricki.—MÉXICO: Michoacán; km marker 81 between Ario de Rosales and La Huacana, 1,602 m, 19°10′59″N, 101°43′42″W (Locality 28; Catalog number not available—TK47890, KX523180; Catalog number not available—TK47897, KX523181*); Las Minas, 3 km SW Tuxpan (Locality 22; Catalog number not available—TK47819, DQ000477; Catalog number not available—TK47807, KX523182*); 13.5 km SW Zitácuaro, UTM 14Q-352122-2140934 (Locality 23; TTU104808, KF201672; TTU104799, KX523183*). Morelos; Cuernavaca, 18°59.142′N, 99°14.130′W, 2210 (Locality 29; BYU20730, KX523184*).

Peromyscus levipes.—MÉXICO: México; 12 km S Acambay (Locality 20; TTU82707, AY322509*); 14.1 km NW Villa del Carbon (Locality 21; TTU90321, KX523178*). Nuevo León; Cola de Caballo (Locality 7; TCWC47956, AF131928*; TTU104866, MW684869); Tlaxcala; 2 km W Teacalco, 2,710 m (Locality 24; TCWC48331, AF131929*).

Peromyscus madrensis.—MÉXICO: Nayarit; Isla María Madre (Locality 14; USNM512599, AF155397*).

Peromyscus schmidlyi.—MÉXICO: Durango; 30 km SW Ojitos (Locality 9; TTU81634, AY322513; Catalog number not available—TK70778, AY322514; TTU115640, AY322515*; TTU81635, AY322520; TTU81638, AY322523); 12 km E Ojitos (Locality 8; TTU81602, AF155405; TTU81610, AY322517; TTU81603, AY322519; TTU81605, AY322522); San Juan de Camarones (Locality 5; TTU81703, AY322524*); 6.1 km W Coyotes, UTM 13-465908E-2634281N (Locality 11; TTU81642, AF155406; TTU81607, AY322516; TTU81611, AY322518; TTU81643, AY322521*; TTU81617, AY370610); Coyotes (Locality 11; ZTNHC: CWK 1997, AF155407; ZTNHC: CWK1993, AF155408). Sonora; 3 km E Yecura, Colegio Yecura (Locality 3; TTU110286, KF201658; TTU110287/TK148275, need GenBank).

Peromyscus simulus.—MÉXICO: Sinaloa; 6.4 km E Concordia, Highway 40 (Locality 12; TCWC45592, AF131927*).

Peromyscus spicilegus.—MÉXICO: Durango; San Juan de Camarones, UTM 13-356961E-2757448N (TTU81640, AY322512*).

Peromyscus stephani.—MÉXICO: Sonora; Isla San Esteban (Locality 1; UMMZ117385, AF155411*).

Peromyscus winkelmanni.—MÉXICO: Michoacán; 6.9 mi WSW Dos Aguas (TCWC45621, AF131930*).

Specimens examined in the karyotypic portion of this study

Peromyscus ensinki.—MÉXICO: Michoacán; 4.8 km S, 3.6 km E Zinapécuaro (Santa Cruz), 14Q-311982-2194265, 2,012 m (Locality 19; TTU104811); 3.5 km S, 4.8 km E Zinapécuaro (Santa Cruz), 14Q-311971-2194257, 2,012 m (Locality 19; TTU110119, TTU110120, TTU110121).

Peromyscus greenbaumi.—MÉXICO: Jalisco; Volcán de Colima (Locality A, see Houseal et al. 1987); Los Reyes (Locality B, see Houseal et al. 1987); Michoacán; 11.8 km WSW Dos Aguas (Locality 26; TCWC45304, TCWC45305); 10.1 km WSW Dos Aguas (Locality 26; TCWC45307, TCWC45309, TCWC45310).

Peromyscus kilpatricki.—MÉXICO: Michoacán; 13.5 km SW Zitácuaro, UTM 14Q-352122-2140934 (Locality 23; TTU104808); 6.6 km S Patzcuaro (Locality C; see Houseal et al. 1987); 14.6 km S Los Azufres (Locality D; TCWC44683, TCWC44684).

Specimens examined in the morphometric portion of this study

Peromyscus beatae.—32 specimens referenced in Schmidly et al. (1988).

Peromyscus boylii.—44 specimens referenced in Schmidly et al. (1988).

Peromyscus carletoni.—54 specimens referenced in Bradley et al. (2014).

Peromyscus greenbaumi.—MÉXICO: Michoacán: 6.3 mi WSW Dos Aguas 7,90s0 ft (circa Localities 26 and 27; TCWC45309–45310); 7.4 mi WSW Dos Aguas 7,900 ft (circa Localities 26 and 27; TCWC45304); 14.0 mi W Aguililla (Dos Aguas) 7,000 ft (circa Localities 26 and 27; UMMZ109536–109537); 2.0 mi W Dos Aguas 7,600 ft (circa Localities 26 and 27; UMMZ109592–109593); Dos Aguas, 7,100 ft (circa Localities 26 and 27; UMMZ109595); 6.3 mi WSW Dos Aguas 8,000 ft (circa Localities 26 and 27; UMMZ109622–109627); 1.0 mi NW Dos Aguas (circa Localities 26 and 27; UMMZ109638); Sierra Varcalosa 1.5 h (by mule) NE Ranch Varcalosa, 8,900 ft (circa Localities 26 and 27; UMMZ102674, UMMZ102675); Michoacán, 8.4 mi N Los Reyes (circa Locality B; TCWC27265).

Peromyscus ensinki.—MÉXICO: Michoacán: 3.5 km S, 4.8 km E Zinapécuaro, 14Q-311982-2194265, 2,012 m (circa Locality 19; TTU110118–110121); 3.6 km S, 4.8 km E Sinecure, 14Q-311971-2194257, 2,012 m (circa Locality 19; TTU114811).

Peromyscus kilpatricki.—MÉXICO: Michoacán: 9.1 mi S Los Azufres (circa Locality D; TCWC284318, TCWC284418); Puerto Garnica (circa Locality D; CWK2032); 2 mi SW Zitácuaro (circa Localities D and 23; KU62864–62865, KU62867–62868,); 13. km SW Zitácuaro (circa Localities D and 23; TTU104799, TTU104808); 3.5 mi N Opopeo (circa Locality C; UMMZ110645); 2.3 mi N Opopeo (circa Locality C; UMMZ110647); 5 mi. S Patzcuaro; 11 km W Quiroga (circa Locality C; UMMZ95714); Uruapan, Cupatitzio National Park (circa Locality E; UMMZ89551, UMMZ89553–89554, UMMZ8955-89567, UMMZ89571-89572, UMMZ89575-89576, UMMZ89578); Parque Nacional, Uruapan (circa Locality E; UMMZ110641–110644); Tzararacua (circa Locality E; UMMZ89581); 5 mi Uruapan (Tzararacua Falls), 5,000 ft (circa Locality E; UMMZ109504, UMMZ109506–109510); Tzararacua Falls, 7 mi S Uruapan 1,400 ft (circa Locality E; UMMZ92126–92134, UMMZ109512–109516); 8.4 mi N Los Reyes (circa Locality B; TCWC27266).

Peromyscus levipes.—85 specimens referenced in Schmidly et al. (1988).

Peromyscus schmidlyi.—9 specimens referenced in Bradley et al. (2004).

Peromyscus simulus.—15 specimens referenced in Schmidly et al. (1988).

Acknowledgments

We thank J. R. Bateman, E. K. Roberts, and S. C. Vrla for reviewing earlier versions of this manuscript. Special thanks to V. J. Swier for assistance in karyotyping some of the specimens. Tissue samples were provided by the Natural Science Research Laboratory, Museum of Texas Tech University, Zadock-Thompson Natural History Collection, University of Vermont, and Monte L. Bean Life Science Museum, Brigham Young University. Collecting permits were provided by SEMARNAT. Special thanks to the Field Methods classes of 2006 and 2008 for assistance in collecting specimens. Portions of this study were supported by a National Science Foundation grant to I. F. Greenbaum, CWK, and DJS, other support for this research was partially provided by a National Institutes of Health Grant (Department of Health and Human Services A141435-01 to RDB).