Exploring species boundaries of the spotted ground
squirrel (Xerospermophilus spilosoma) complex

Sergio Albino-Miranda1,2 , Juan M. Uriostegui-Velarde3 , Pablo F. Colunga-Salas1 , Yessica Rico4 ,
Alejandro A. Castro-Luna1 , Alberto González-Romero2 , and Jorge Galindo-González1* .

1Instituto de Biotecnología y Ecología Aplicada (INBIOTECA), Universidad Veracruzana. Campus para la Cultura, las Artes y el Deporte, Av. De las Culturas Veracruzana No. 101, CP. 91094, Zona Universitaria, Xalapa, Veracruz, México. pcolunga@uv.mx, (PFC-S); alcastro@uv.mx, (AAC-L).

2Red de Biología y Conservación de Vertebrados, Instituto de Ecología, A.C. (INECOL). Carretera Antigua a Coatepec No. 351, CP. 91073, El Haya, Xalapa, Veracruz, México. sergio.albino@inecol.mx, (SA-M); alberto.gonzalez@inecol.mx, (AG-R).

3Facultad de Ciencias Biológicas, Universidad Autónoma del Estado de Morelos. Avenida Universidad No. 1001, CP: 62209, Chamilpa, Cuernavaca, Morelos. juan.uriosteguive@docentes.uaem.edu.mx, (JMU-V).

4Red de Diversidad Biológica del Occidente Mexicano, Centro Regional del Bajío Pátzcuaro, Michoacán, Instituto de Ecología, A.C. (INECOL). Av. Lázaro Cárdenas No. 253 CP. 61600, Centro, Pátzcuaro, Michoacán. yessica.rico@inecol.mx, (YR).

*Corresponding author: jgalindo@uv.mx

Xerospermophilus spilosoma exhibits notable geographic and morphological variation, prompting a debate over its taxonomic status. Currently, it is unclear whether it represents a single highly variable species or a complex of cryptic species that includes X. perotensis. Although the latter has a larger body size and a distinctive dorsal pattern, current genetic analyses do not support its recognition as a separate species. This study aimed to delimit potential evolutionary units within the X. spilosoma complex. Twenty-four sequences of the mitochondrial cytochrome b gene were analyzed using 690-bp fragments from X. spilosoma and X. perotensis specimens collected from eight locations. Phylogenetic and divergence time inferences were estimated using Maximum Likelihood and Bayesian Inference, along with analyses of genetic distances and haplotype networks. Three species delimitation methods (ABGD, PTP, and GMYC) were applied, and the ecological uniqueness and areas of overlap within and between species of the X. spilosoma complex were assessed. Four lineages comprising 22 unique haplotypes were identified, with interpopulation genetic distances ranging from 3 % to 6 %. Species delimitation methods suggested between one and four potential species. Meanwhile, the comparison of ecological niches revealed limited overlap. Genetic and environmental evidence indicate that X. spilosoma comprises at least three evolutionarily independent lineages. The complex originated in the Miocene, more than 5 million years ago, with divergence events concentrated between 3.5 and 1.5 million years ago, in accordance with geographical barriers such as Río Grande and the Nazas River and the Trans-Mexican Volcanic Belt. These results highlight the need to conserve these populations as independent evolutionary units, particularly the Perote population, given its isolation and ecological and genetic uniqueness.

Keywords: Chihuahuan Desert; cytochrome b; Río Grande; Taxonomy.

Xerospermophilus spilosoma presenta una notable variación geográfica y morfológica, lo que ha generado debate sobre su estatus taxonómico. Actualmente se discute si representa una única especie o un complejo de especies crípticas que incluiría a X. perotensis. Aunque esta última presenta mayor tamaño corporal y un patrón dorsal distintivo, los análisis genéticos actuales no respaldan su reconocimiento como especie válida. El objetivo de este trabajo fue delimitar las posibles unidades evolutivas dentro del complejo X. spilosoma. Se analizaron 24 secuencias del gen mitocondrial citocromo b, usando fragmentos de 690 pb de X. spilosoma y X. perotensis, de ocho localidades distintas. Las inferencias filogenéticas y de tiempo de divergencia se estimaron mediante Máxima Verosimilitud e Inferencia Bayesiana, junto con como análisis de distancias genéticas y redes de haplotipos. Se aplicaron tres métodos de delimitación de especies (ABGD, PTP y GMYC), además de evaluarse el nicho ecológico por especie y sobreposicionamiento entre especies del complejo X. spilosoma. Se identificaron cuatro linajes con 22 haplotipos únicos con distancias genéticas interpoblacionales entre el 3 y el 6%. Los métodos de delimitación de especies sugieren entre una y cuatro especies potenciales. Por su parte, la comparación de los nichos ecológicos mostró un bajo solapamiento de áreas. La evidencia genética y ambiental obtenida sugiere que X. spilosoma corresponde a un complejo de al menos tres linajes evolutivamente independientes. El complejo se originó en el Mioceno, hace más de 5 Ma, con eventos de divergencia concentrados entre 3.5 y 1.5 Ma, en concordancia con barreras geográficas como los ríos Grande y Nazas y la Faja Volcánica Transmexicana. Estos resultados, resaltan la necesidad de conservar estas poblaciones como unidades evolutivas independientes, especialmente la de Perote, por su aislamiento y singularidad ecológica y genética.

Palabras clave: Citocromo b; Desierto Chihuahuense; Río Grande; Taxonomía.

Establishing the boundaries between species is a difficult task, especially in taxa for which speciation processes have not resulted in an evident morphological differentiation (Goldstein and De Salle 2011; Fišer et al. 2018). This can be explained by phenomena such as convergence (Losos 2008, 2011), stabilizing selection (Gould 2002; Hansen and Houle 2004), or a recent speciation process (Gittenberger 1991; Rundell and Price 2009). A major problem in taxonomy is that species boundaries vary widely depending on the species concept employed (De Queiroz 2007). Therefore, a unified concept based on their evolutionary origin has been proposed, defining them as evolutionarily independent lineages. These lineages are identified by evaluating secondary characteristics, e.g., ecological, morphological, and genetic traits, that reflect some degree of evolutionary independence (Simpson 1961; Wiley 1978; Bock 2004; Hey 2006; De Queiroz 2007).

With more than 2600 species, the order Rodentia is the most diverse group of mammals worldwide (D’Elía et al. 2019). This diversity is partly due to their high evolutionary rates and rapid radiation processes, often associated with evolutionary convergences, which has produced taxonomic complexities that still require resolution (Triant and De Woody 2006; Fabre et al. 2012; Burgin et al. 2018). Within this order, squirrels of the family Sciuridae stand out for their diversity, with about 300 species distributed in a wide range of ecosystems, from deserts to tropical forests, and from sea level to more than 4500 meters above sea level. This family exhibits great morphological variation, including arboreal, terrestrial, and gliding forms, with solitary or social patterns (Koprowski et al. 2016; Rocha et al. 2016). In addition, their biological characteristics, such as low dispersal capacity, short life cycles, and varied reproductive strategies, allow for a detailed analysis of the differentiation and evolutionary isolation processes, making squirrels ideal models for genetic and evolutionary studies (Rocha et al. 2016; Flores-Manzanero and Vázquez-Domínguez ٢٠١٩; Waterman et al. 2021).

The current distribution of most squirrel species in North America is explained by allopatric speciation processes, where geographic barriers such as mountain systems and water bodies have limited gene flow (Harrison et al. 2003; Ge et al. 2014; Zelditch et al. 2015). For example, the Snake River, in the northwestern United States, has influenced the diversification of small-eared ground squirrel species of the genus Urocitellus (McLean et al. 2025). In this context, several current species have been isolated in valleys due to geological events that occurred during the Quaternary period, such as fluctuating glaciations or volcanic activity (Harrison et al. 2003; Van Tuinen et al. 2008; Menéndez et al. 2021). A peculiar case is that of the Mojave ground squirrel, Xerospermophilus mohavensis, whose possible origin involved allopatric speciation in a small isolated refuge within the Mojave Desert, delimited by the Sierra Nevada (Bell et al. 2010).

However, geographic isolation does not always imply sufficient ecological or genetic divergence to justify the separation into distinct species (Barton 2020; Kulmuni et al. 2020). An example is the population of Mearns squirrels (Tamiasciurus mearnsi) in Baja California, for which genetic studies indicate a relatively recent separation and a limited ecological differentiation despite being geographically isolated from Tamiasciurus douglasii, leading to the question of whether T. mearnsi should really be considered a distinct species (Arbogast et al. 2001; Pecnerová and Martinova 2012; Hope et al. 2016). Therefore, examining the niches of species in the environmental space and projecting them in the geographical space constitutes an additional line of evidence to evaluate whether populations, in addition to being physically separated, have occupied and exploited different environments, which could indicate processes of ecological differentiation that support or facilitate the delimitation of species (Peterson et al. 2000; Kozak and Wiens 2006; Pahad et al. 2019). Similarly, niche modeling studies contribute to understanding the relationship between genetic variation patterns and the environmental parameters that limit the distribution of a species (Ashrafzadeh et al. 2018; Calixto-Pérez et al. 2018; Luna-Aranguré and Vázquez-Domínguez ٢٠٢٠).

A taxonomic review of the genus Spermophilus based on morphology and sequences of the cytochrome b (cytb) mitochondrial gene revealed its paraphyletic condition, which led to the reassignment of several species to new genera. Thus, Xerospermophilus was recognized as a genus that includes X. mohavensis, X. tereticaudus, X. spilosoma, and X. perotensis (Helgen et al. 2009). However, the phylogenetic relationship between X. spilosoma and X. perotensis remains a matter of debate, particularly with regard to their recognition as separate species (Harrison et al. 2003; Helgen et al. 2009; Fernández 2012).

Xerospermophilus spilosoma is the most divergent species that is distributed over a wide geographic range, from Wyoming, South Dakota, and Nebraska to central Mexico (Lacher et al. 2016; Álvarez-Castañeda 2024). This wide distribution has favored a high morphological diversity within the species, including variation in coat coloration from cinnamon to dark smoky-gray shades, through various shades of brown. In addition, it has a series of spots on the back and flanks whose shape and intensity vary depending on the environment, probably in response to local ecological conditions (Cothran et al. 1977; Álvarez-Castañeda 2024). As a result of this morphological variation, 13 nominal subspecies have been described (Helgen et al. 2009; Mammal Diversity Database 2025).

The taxonomic classification of X. spilosoma is currently a subject of debate. It is not clear whether it is a single species with high population variability that includes several subspecies (Helgen et al. 2009; Fernández 2012), including X. perotensis, or if the possible paraphilia observed between the subspecies of X. spilosoma and X. perotensis indicates that it could be a complex of species (Harrison et al. 2003).

Although X. perotensis has a larger body size, differences in coloration with a distinctive pattern of dorsal spots and distinct vocalizations, the genetic analyses available to date, based on nuclear (GHR and IRBP) and mitochondrial (Cytb and 12S rRNA) genes, suggest that genetic divergence would not be sufficient to support its recognition as a separate species (Fernández 2012). However, it is worth noting that these conclusions were based on an limited number of specimens, namely 3 individuals of X. spilosoma and 4 of X. perotensis. Therefore, the objective of this study was to evaluate the differentiation of populations within the X. spilosoma complex and their possible status as species using four distinct lines of evidence: (1) phylogenetically delineating the evolutionary units; (2) estimating genetic distances between populations to quantify their divergence; (3) applying species delimitation methods that allow the identification of separate evolutionary lineages; and (4) comparing the degree of ecological differentiation between populations.

Materials and methods

Sample collection. A total of 23 tissue samples were processed, from specimens collected in the field (n = 12) and from specimens obtained from different national and international scientific collections (n = 11). The field collection consisted of ectomization of the third phalanx from X. perotensis and X. spilosoma specimens captured in Perote and Mapimí, respectively. Sample collection followed the recommendations of Romero-Almaraz et al. (2007) and the guidelines of the Committee on Animal Care and Use of the American Society of Mammalogists (Sikes et al. 2016). The corresponding collection permits were issued by the Board of the Environment and Natural Resources (SPARN/DGVS/04074/23 and SPARN/DGVS/08359/23). The tissues were fixed separately in 96% ethanol and stored at -70 °C until processing.

Genomic DNA (gDNA) was extracted from the tissues using the DNeasy Blood and Tissue Kit (QIAGEN, USA) according to the manufacturer’s instructions. Subsequently, the concentration of the genetic material was quantified in a Nanodrop spectrophotometer (model DS-11, Denovix). To perform the amplification of cytb fragments by Polymerase Chain Reaction (PCR), a new pair of oligonucleotides was designed in PerlPrimer v1.1.21: Xeros_Fw (5’YSAYTTACMYGCACCYTCC-3’) and Xeros_Rv (5’GGRTATWCAACRGGTTGYCMTC 3’), which amplifies a 981-bp fragment. This primer pair was evaluated in silico using SnapGene v.8.1 to verify its specificity, hybridization efficiency, and absence of secondary structures.

The PCR reactions were run in a final volume of 25 μL per sample, containing 2 μL of gDNA (with concentrations between 20 and 200 ng/μL), 14.25 μL of nuclease-free H2O, 6.25 μL of Master mix DreamTaq (QIAGEN, USA), and 2 μL of each primer at 10 μM. One positive control and one negative control were included in each run to ensure procedural reliability and to rule out contamination. The amplification protocol consisted of an initial denaturation step at 94 °C for ٥ minutes, followed by 40 denaturation cycles at 94 °C for 30 seconds, alignment at 56 °C for 30 seconds, extension at 72 °C for 70 seconds, and a final elongation step at 72 °C for four minutes. The PCR tests were performed on a BioRad T100 thermal cycler. PCR products were visualized on 2 % agarose gels, stained with red gel, and examined under UV light using a transilluminator. The PCR products were then purified using the Qiaquick Gel Extraction Kit (QIAGEN, USA) according to the manufacturer’s instructions. The purified products were sent for bidirectional sequencing to Macrogen, South Korea.

In total, 24 Xerospermophilus samples were considered for subsequent analyses. Of these, 12 were generated in this study (7 for X. spilosoma and 5 for X. perotensis), and 12 were recovered from GenBank. The sequences were edited manually, then aligned using the MUSCLE algorithm, and visually inspected using Geneious Prime and MEGA X (Kumar et al. 2018). The final alignment included 14 sequences for X. spilosoma, 10 for X. perotensis (Appendix 1), and one for Ictidomys mexicanus as an outgroup (Harrison et al. 2003; Guevara-Chumacero et al. 2006; Fernández 2012).

Phylogenetic analyses and genetic diversity. From the global alignment, the best substitution model was selected through an exhaustive search in JModelTest v2.1.10 (Darriba et al. 2012). Phylogenetic inference was performed under the Maximum Likelihood (ML) criterion using the previously selected model (GTR + G), with estimated branch support across 10 000 ultrafast bootstrap replicates (Hoang et al. 2018) in IQ-TREE v1.6.12 (Nguyen et al. 2015). Bayesian Inference (BI) was performed in MrBayes v.3.2.7 (Ronquist et al. 2012) with the MCMC algorithm and the previously calculated substitution model (GTR + G + I); to this end, two separate runs were performed with three hot chains and one cold chain with 10 million generations and a 25 % burn-in. Chain convergence and good sampling (ESS > 200) were assessed in TRACER v1.7 (Rambaut et al. 2018). The consensus trees generated in each inference were visualized and edited using FigTree v1.4.4. (Rambaut 2010).

Subsequently, a haplotype file was generated in DnaSP v6 (Rozas et al. 2017), from which the number of haplotypes was identified, and haplotype diversity (Hd) and nucleotide diversity (Pi) were calculated. With this information, a haplotype network was constructed using the TCS criterion in PopArt v1.7 (Leigh et al. 2015). Additionally, an analysis of intra- and interclade genetic distances was performed using uncorrected pairwise sequence distances (p-distances) using the ape library (Paradis et al. 2004). Based on these data, a heat map was drawn up with the adegenet and split libraries (Jombart 2008; Ezard et al. 2009), both in R v4.4.1 (R Core Team 2020). Genetic distances within the genus Xerospermophilus were evaluated using 15 sequences of X. tereticaudus and 10 sequences of X. mohavensis (Harrison et al. 2003; Bell et al. 2010; Fernández 2012).

Divergence times. With the alignment constructed as described above, divergence times were inferred in BEAST v2.6 (Bouckaert et al. 2019), considering a non-correlated lognormal relaxed molecular clock and replacing I. mexicanus for Cynomys ludovicianus (Appendix 1) as the outgroup (Castellanos-Morales et al. 2014). GTR substitution models with the gamma distribution and the Bayesian Coalescent Skyline plot model were established.

A secondary calibration based on four nodes was performed. The first node was calibrated based on the fossil record of Cynomys rafinesque, with a maximum age of 1.8 Ma (Ge et al. 2019), representing the divergence between Cynomys and the X. spilosoma complex. The second node corresponds to the separation between the Kansas, Texas, and New Mexico populations, in the USA, and those of Arizona and Mexico, attributed to the formation of the Rio Grande, whose average consolidation is estimated to have occurred 2.6 Ma ago (Morgan and Golombek 1984; Repasch 2017). The third node reflects the divergence between the populations of Arizona, USA, and Mapimí, Mexico, in contrast to the populations of the Central Plateau of Mexico, associated with the Nazas River, whose origin dates back to the Pliocene and with its current configuration estimated to have occurred 1.5 Ma ago (Petersen 1976; Hafner and Riddle 2011). Finally, the fourth node represents the isolation of the Perote population, in Mexico, attributable to the culmination of the formation of the Trans-Mexican Volcanic Belt, completed approximately 1.6 Ma ago (Ferrari et al. 1999, Ferrari 2000; Gámez et al. 2017).

All nodes were assigned a lognormal distribution as a prior, as this model reflects the observation that evolutionary divergence typically precedes both the first fossil record and the complete establishment of a geographic event. Under this perspective, speciation does not occur at the same time that a fossil appears or a barrier is formed; instead, these elements represent only minimum age limits (Ho 2007). MCMC analyses were performed with two independent runs of 10 million generations each, sampling every 1000 generations. The convergence, stability, and adequate sampling of the results were evaluated with Tracer v1.7.1 (Rambaut et al. 2018), applying a 25 % burn-in. The final phylogenetic tree that included the divergence time intervals was generated using TreeAnnotator v2.6.2 (Bouckaert et al. 2019). The graphs were generated using the deeptime (Gearty 2025) and ggplot2 (Wickham 2016) packages in R.

Species delimitation. Boundaries between species within the X. spilosoma complex were evaluated using three delimitation methods: the Automatic Barcode Gap Discovery (ABGD), which is based on the identification of a barcode corresponding to the natural discontinuity between intraspecific (minor) and interspecific (major) genetic variability. ABGD groups sequences according to this discontinuity, without the need for an a priori hypothesis about the number of species (Puillandre et al. 2012). The Poisson Tree Processes (PTP) method, based on a maximum likelihood model, assesses differences in substitution rates between clades to detect clusters that are consistent with putative species (Zhang et al. 2013), while the Generalized Mixed Yule Coalescent (GMYC) method uses an ultrametric tree to analyze branch lengths and determine whether they correspond to intraspecific or interspecific processes (Fujisawa and Barraclough 2013).

For ABGD, the previously calculated genetic distance matrix was used. The configuration of the parameters was (i) a minimum intraspecific distance (Pmin) of 0.001 and a maximum intraspecific distance (Pmax) ranging from 0.02 to 0.1; (ii) a Barcode Gap Width of 1.5; and (iii) the Jukes-Cantor model (JC69); for the analysis, we used the command line in the ABGD program (Puillandre et al. 2012). For PTP, we used the ML tree generated with 100 000 MCMC generations and a 10 % burn-in on the PTP web server (http://species.h-its.org/). For GMYC, we used an ultrametric guide tree of the genus Xerospermophilus generated in MEGA X (Kumar et al. 2018). A multiple method was applied using a λ = 5 to fit the molecular-clock model and evaluated confidence intervals (0 - 2); the model was performed using the splits package in R (Ezard et al. 2009).

Niche modeling. In addition to the genetically analyzed samples, we included 404 records of the presence of X. spilosoma and X. perotensis from GBIF (2024). These records were organized according to the identification of terrestrial ecoregions proposed by Olson et al. (2001) corresponding to each clade defined in the phylogenetic analyses and the delimitation of potential geographical barriers between records: the Rio Grande between clades A and B, the Nazas River between clades B and C, and a minimum convex polygon for clade D, subsequently intersected with the aforementioned ecoregions. A 6 km buffer was established for each polygon delimited by these ecoregions and barriers. This distance was set given the low mobility reported for these squirrels, whose mean dispersal distance is approximately 1.5 km (maximum 2.8 km) (Montero-Bagatella and González-Romero 2014). For clades A and D, the delimitation was obtained directly with this procedure. For clade B, a northern limit was established considering a 30 km buffer with respect to the clade A polygon and applying a 21 km reduction to the south. For clade C, a 21 km reduction to the north was applied, maintaining the previous delimitation to the south, which delineated a 30 km strip between the clades A, B, and C polygons. Based on these polygons, records were classified and grouped for modeling by clades, excluding those located within the 30 km strips or outside the polygons; this procedure was carried out in ArcMap v.10.5.

The ecological niche models were constructed using the 19 bioclimatic variables available in the WorldClim platform (WorldClim 2024), which are in raster format with a 30 arcsec (1 km²) resolution. The models were generated using the Maximum Entropy algorithm, implemented in MAXENT v3.4.1 (Phillips et al. 2017). First, a general model was constructed for each clade to identify the most representative variables. To this end, we used the contribution and permutation percentage table, as well as the Jackknife test, both available in the MAXENT outputs, to assess the relative importance of the explanatory variables. Based on the variables selected by clade, a consensus was performed, selecting those relevant to the four clades.

Subsequently, a Pearson correlation analysis was performed on the previously selected variables. Highly correlated variables (r > 0.85) were excluded to avoid collinearity and reduce redundant information (Segurado et al. 2006). We selected seven correlated variables with a simpler environmental interpretation that do not combine humidity and temperature data: BIO1 (mean annual temperature, °C), BIO2 (mean diurnal temperature range), BIO4 (temperature seasonality, %), BIO10 (mean temperature of the warmest quarter, °C), BIO11 (mean temperature of the coldest quarter, °C), BIO15 (seasonality of precipitation, variability index), and BIO16 (precipitation of the wettest quarter, mm).

Once the clade-based groups were defined, we performed a 1 km spatial filtering of the final records (clade A, 151; clade B, 72; clade C, 37; clade D, 37) using the Wallace package (Kass et al. 2018). Then, the models for each group were calibrated using the accessibility or mobility area with 20 000 randomly selected background pixels. This selection represents a hypothesis about the area to which the species has or has had access to disperse (Barve et al. 2011). Subsequently, a spatial partitioning of occurrence data was performed to train and validate the models. The method implemented was the chessboard, considering an aggregation factor of 2 (Muscarella et al. 2014). Models with linear, quadratic, and combined functions were fitted using regulation multipliers between 0.5 and 3 (0.5 intervals), selecting the best models according to the AIC (Burnham and Anderson 2002). From these, binary (presence/absence) distribution maps were generated according to the biogeographic provinces of Olson et al. (2001).

Additionally, niches were characterized in a multivariate space using a Principal Component Analysis (PCA), performing the ordination through a correlation matrix and using Mahalanobis distance approximations (Broennimann et al. 2012). The degree of overlap or divergence between the ecological niches of the candidate species was evaluated using similarity and equivalence statistical tests implemented in the ecospat v4.2.1 package (Di Cola et al. 2017). In both cases, the overlap between niches was quantified using Schoener’s-D index, whose values range from 0 (completely discordant niches) to 1 (identical niches) (Schoener 1970; Warren et al. 2008).

The niche similarity analysis assesses whether the niches of two species are more similar than expected by chance, accounting for the environmental context in which they occur. In this test, a P-value < 0.05 suggests that niches are more similar than expected by chance (Broenniman et al. 2012). In contrast, in the niche equivalence test, the observed D-value is compared to a null distribution constructed from random reassignments of occurrences of each species pair (Brown and Carnaval 2019). In this test, a P-value < 0.05 indicates that the niches are more different than expected by chance. For each test, 100 replicates were performed to generate a null distribution of overlap scores, which was compared to the observed values.

Results

The final alignment comprised 690 bp, with 146 variable sites, 81 of which were parsimony-informative. The ML and BI topologies were consistent, recovering the same phylogenetic relationships and with support values greater than 0.85 / 85 %, respectively. Four clades were identified (Figure 1a): clade A includes sequences from the USA populations (Kansas, Texas, and New Mexico); clade B, from Arizona, USA, and Durango, Mexico; clade C only includes sequences from San Luis Potosí, Mexico; finally, clade D includes sequences from the Perote valley in Puebla and Veracruz, Mexico.

A total of 22 unique haplotypes were recovered, with a haplotype diversity (Hd) of 0.989 and a nucleotide diversity (Pi) of 0.04360. The haplotype network showed four clusters corresponding to the clades observed in the phylogenetic relationships (Figure 1c). These clusters were separated by different numbers of mutational steps (7-14), indicating clear genetic differences between them. Likewise, a high number of mutational steps were observed between the Arizona and Mapimí populations, as well as between the New Mexico and Texas–Kansas populations. In addition, several hypothetical haplotypes were identified that connect the various clusters. On the other hand, the genetic distances between I. mexicanus and the remaining sequences analyzed were low (10–11 %). In X. spilosoma, intrapopulation genetic variability ranged from 1 % to 3 %, whereas interpopulation variability ranged from 3 % to 6 % (Figure 2).

The estimated divergence times suggest that the X. spi-losoma complex originated toward the end of the Miocene, just over 5 Ma ago (Figure 3). Most divergence events are concentrated between 3.5 and 1.5 Ma, with a marked increase in the divergence process during the Pliocene.

Of the eight partitions generated in the ABGD analysis, the three-species option yielded the highest statistical significance (P = 0.007; Figure 1b). On the other hand, the PTP analysis identified between 2 and 6 potential species; however, only the four species with support values close to 0.5 were considered. Finally, the GMYC clustering analysis identified only one evolutionary entity, with a likelihood ratio of 19.12 (P = 7.04e-05).

Ecological niche analyses using PCA indicate consistent ecological segregation across most clades, except for clades A and B, which show a smaller Mahalanobis distance between their environmental centroids (Figure 4). The overlap between clades is low to moderate, suggesting possible ecological differentiation (Figure 5). In contrast, clade A and clade B showed a partial differentiation with some overlap (D = 0.23). The comparisons of clades B–C and B–D showed the least separation, with overlap approaching 1 in both comparisons (D = 0.92). The analysis indicated no evidence of similarity greater than expected by chance (P > 0.05), with the exception of clades B vs. D (P = 0.03; Table 1). On the other hand, in the equivalence test, none of the comparisons showed significant differences (P < 0.05), indicating that the niche models generated for each clade pair did not differ in their niches within the shared environmental space, with D-index values ranging from 0 to 0.14 (Table 1).

Discussion

The results suggest that Xerospermophilus spilosoma comprises a complex that hosts four evolutionarily independent lineages, a hypothesis supported by clear phylogenetic separation, consistent haplotype distributions across environmentally segregated populations, and differentiation of ecological niches among these lineages.

The estimated divergence times differ from those proposed by Fernández (2012), who suggested that the X. spilosoma complex is approximately 2.69 Ma old and that the last common ancestor of the San Luis Potosí and Perote lineages existed 0.74 Ma ago. In contrast, our results indicate that the X. spilosoma complex emerged toward the end of the Miocene, just over 5 Ma ago. Most divergence events are concentrated between 3.5 and 1.5 Ma and are consistent with the estimated ages of the geographical barriers that influenced their diversification.

Clade A is the first to diverge, possibly favored by the geographical barrier represented by the Rio Grande (Figure 6). This river, with a length of approximately 3,050 km, flows through a large portion of the southwestern United States, emptying into the Gulf of Mexico (Kelley 1952). It was formed during the Late Miocene and Pliocene between 6.9 and 2.5 Ma (Morgan and Golombek 1984; Gamez et al. 2017), consistent with the divergence time estimated in this work. The importance of this water body as a driver of speciation processes has been documented for other rodents, such as Chaetodipus and Eutamias (Sullivan 1985; Neiswenter et al. 2019).

Clade B is distributed in the Trans-Pecos and Bolsón de Mapimí, while clade C is represented by the population inhabiting the Central Plateau in San Luis Potosí, Mexico. These three zones are subregions of the great Chihuahuan Desert, considered among the most diverse arid regions worldwide (Dinerstein et al. 2000; Edwards et al. 2001). The geological and climatic events of the Miocene and Pliocene, along with the cyclical climatic changes of the Pleistocene, had a profound effect on the genetic structure and distribution of various species that live in this desert (Raymo and Ruddiman 1992; Hafner and Riddle 2005; Loera et al. 2017; Scheinvar et al. 2020). These events possibly favored divergence between clades, which is supported by the divergence times estimated in this study, similar to the reports regarding the diversification between Cynomys mexicanus and C. ludovicianus (Castellanos-Morales et al. 2016), as well as about the genetic patterns of Perognathus and Chaetodipus mice (Riddle et al. 2000; Neiswenter and Riddle 2010) and the differentiated distributions of grasshopper mice of the genus Onychomys (Riddle 1995).

On the other hand, the divergence between clades B and C could be related to the Nazas River, whose formation is estimated to have occurred between mid- and late Pliocene (Petersen 1976; Hafner and Riddle 2011), a period estimated in the present work. This river has been identified as a possible biogeographic barrier for several species, including the rats Neotoma albigula and N. leucodon (Edwards et al. 2001), the cactus mouse Peromyscus eremicus (Riddle et al. 2000), gophers of the genus Cratogeomys (Hafner et al. 2008), and the black-tailed hare Lepus californicus (Lorenzo et al. 2021).

On the other hand, clade D is restricted to the Eastern Basin, a high-mountain region, formed mostly during the Pleistocene (<1.6 Ma) as part of the uplift of the Trans-Mexican Volcanic Belt (Ferrari et al. 1999). The emergence of this mountain range could have favored the divergence of this clade, as documented for other rodents, such as Dipodomys phillipsi and Peromyscus bullatus (Peterson et al. 2000; González-Ruíz et al. 2005; Sánchez-Cordero et al. 2005; Fernández et al. 2012).

The interspecific genetic distances observed within the X. spilosoma complex were slightly below the threshold proposed by Baker and Bradley (2006) to recognize entities as distinct species. However, these distances are similar to those accepted for other closely related genera, such as Spermophilus, when supported by additional evidence, such as genetic structure and multilocus analysis (Simonov et al. 2024). On the other hand, these divergence levels are consistent with evolutionarily independent lineages, as proposed for squirrel species of the genus Tamias, which supports their possible recognition as distinct species (Ge et al. 2014). The genetic diversity patterns found in X. spilosoma in the present study reinforce this interpretation. A high Hd = 0.99 and a considerable Pi = 0.04 were observed, which are within the upper range reported for other rodents, such as Neotoma mexicana (Hd = 0.97; Pi 0.03–0.05) and Reithrodontomys chrysopsis (Hd = 0.99; Pi = 0.03) (Hernández-Canchola et al. 2021; León-Tapia et al. 2023).

On the other hand, the different species delimitation methods yielded inconsistent results, as they identified different numbers of evolutionary entities. Such discrepancies between methods are common, as each approach is based on different assumptions and models (Feijó et al. 2019; Martínez‐Borrego et al. 2023). The partial concordance between ABGD and PTP suggests some degree of genetic structure within X. spilosoma, although the low resolution of GMYC could indicate recent diversification or insufficient mitochondrial differentiation, which warrants further analysis.

In addition, correspondence between genetic structure and ecological niche modeling was observed. Clades A, C, and D exhibited marked ecological differentiation, with minimal or no niche overlap among them. In contrast, we found evidence of partial or complete niche overlap between clade B and all other clades. This overlap can be explained by the niche breadth of clade B, associated with the marked climatic heterogeneity of the Mapimí region, which offers a wide range of environmental conditions (Van Devender and Burgess 1985; García-Arévalo and Nocedal 2008). Although the overall pattern is a positive relationship between niche breadth and geographic range (Morin and Lechowicz 2013; Slatyer et al. 2013; Moulatlet et al. 2025), niche breadth does not necessarily correspond to a larger potential distribution area, since the conditions that make up that niche may be unequally represented in geography (Dallas and Ten Caten 2025). For example, clade A has a potential distribution area of 637 867 km², more than twice the estimated area for clade B (309 468 km²), with a smaller environmental niche breadth; this highlights the importance of the environmental heterogeneity of the Mapimí region and is related to the distribution of clade B (García-Arévalo and Nocedal 2008).

Statistically, the ecological niche of clade B showed similarity only with the niche of clade D. Niche similarity between clades B and D could be explained by a non-typical niche phylogenetic conservatism, in which species maintain environmental similarities regardless of their genetic distances (Wiens and Graham 2005: Losos 2008). It is also possible that both clades are preserving a niche close to a midpoint within the environmental space occupied by the clade set (Wiens et al. 2010; Peixoto et al. 2017). However, the geographical distance between populations, genetic evidence, and the limited dispersal capacity of squirrels suggest that clades B and D are distinct groups despite their ecological proximity. The equivalence test was not significant, implying that the compared niches cannot be considered identical. This is consistent with the fact that the equivalence test is more conservative than the similarity test, as it assesses only whether two niches are indistinguishable from actual locations, without incorporating the surrounding environmental space (Brown and Carnaval 2019).

Taken together, the lines of evidence discussed here highlight the importance of considering both evolutionary history and local ecological conditions in order to understand diversification processes within the X. spilosoma complex. The evidence from the present study suggests that this complex can be divided into four lineages. Since species delimitation should not be based solely on absolute genetic distances, it is essential to integrate other factors such as monophyly, phylogenetic structure, and the gradual nature of allopatric speciation, the pace of which may vary according to the ecological, genetic, and geographical conditions involved (Carstens et al. 2013). In addition, it is important to recognize that genetic divergence can advance even in the absence of differentiated ecological selective pressure (Wiens 2004; Nosil 2012).

In the case of clade D (X. perotensis), the data support its recognition as an evolutionarily distinct lineage characterized by geographic isolation, a unique genetic structure, monophyly, and ecological differentiation from the other clades. These elements support its recognition as a sister species of X. spilosoma under the unified species concept framework. However, inconsistencies between the delimitation methods and the molecular evidence, limited to a mitochondrial gene fragment, preclude a definitive taxonomic conclusion. An alternative interpretation is that X. perotensis is at an advanced stage of the allopatric speciation continuum.

Although only a partial region (690 bp) of the cytochrome b gene was analyzed, our results provide solid evidence on genetic structuring and ecological differentiation in squirrels of the genus Xerospermophyllus. A more complete and robust reconstruction of the evolutionary history of this taxonomic complex warrants the incorporation of additional mitochondrial and nuclear markers, along with other types of evidence (e.g., morphological, behavioral, or genomic). It will also be essential to increase the sample size and include additional populations, such as those inhabiting the Great Tamaulipas Desert, which could not be included in this study due to insufficient data. Finally, a major limitation in the reconstruction of ecological niche models is the lack of biotic data, such as ecological interactions, resource availability, or predator pressure, which restricts the power of models to more accurately discriminate between ecological niches of different clades (Peng et al. 2025).

While the taxonomy of this squirrel complex is being resolved, it is important to recognize that the identification of genetically differentiated clades and distinct ecological niches underscores the need to conserve each genetic lineage as a significant evolutionary unit. This is particularly relevant for the Perote population, the southernmost and most isolated squirrel population, which, due to its distinctive characteristics, is essential for preserving the genetic and ecological diversity of this group.

Acknowledgements

We thank J. Ramos-Luna and E. López-Tello Mera for their support during fieldwork, and C. Flores for supporting the primer design. We also thank the curators for facilitating access to tissue samples for genetic analysis, including T. Lee from the Natural History Collection of the Christian University of Abilene (ACUNHC), F. Cervantes and Y. Hortelano from the National Collection of Mammals, Universidad Nacional Autónoma de Mexico (CNMA), C. López and D. García from the Mammal Collection of the CIIDIR. Secretariat of Science, Humanities, Technology and Innovation (SECIHTI), for the doctoral scholarship number 4001596 awarded to Sergio Albino Miranda.

Literature cited

Álvarez-Castañeda ST. 2024. Mammals of North America-Systematics and Taxonomy. Vol. 2. Cham (CHE): Springer. https://doi.org/10.1007/978-3-031-50825-7

Arbogast BS, Browne RA, and Weigl PD. 2001. Evolutionary genetics and Pleistocene biogeography of North American tree squirrels (Tamiasciurus). Journal of Mammalogy 82:302–319. https://doi.org/10.1644/1545-1542(2001)082<0302:EGAPBO>2.0.CO;2

Ashrafzadeh MR, Khosravi R, Ahmadi M, and Kaboli M. 2018. Landscape heterogeneity and ecological niche isolation shape the distribution of spatial genetic variation in Iranian brown bears, Ursus arctos (Carnivora: Ursidae). Mammalian Biology 93:64–75. https://doi.org/10.1016/j.mambio.2018.08.007

Baker RJ, and Bradley RD. 2006. Speciation in mammals and the genetic species concept. Journal of Mammalogy 87:643–662. https://doi.org/10.1644/06-MAMM-F-038R2.1

Barton NH. 2020. On the completion of speciation. Philosophical Transactions of the Royal Society B 375:20190530. https://doi.org/10.1098/rstb.2019.0530

Barve N, Barve V, Jiménez-Valverde A, Lira-Noriega A, Maher SP, Peterson AT, et al. 2011. The crucial role of the accessible area in ecological niche modeling and species distribution modeling. Ecological Modelling 222:1810–1819. https://doi.org/10.1016/j.ecolmodel.2011.02.011

Bell KC, Hafner DJ, Leitner P, and Matocq MD. 2010. Phylogeography of the ground squirrel subgenus Xerospermophilus and assembly of the Mojave Desert biota. Journal of Biogeography 37:363–378. https://doi.org/10.1111/j.1365-2699.2009.02202.x

Bock WJ. 2004. Species: the concept, category and taxon. Journal of Zoological Systematics and Evolutionary Research 42:178–190. https://doi.org/10.1111/j.1439-0469.2004.00276.x

Bouckaert R, Vaughan TG, Barido-Sottani J, Duchêne S, Fourment M, Gavryushkina A, et al. 2019. BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis. PLoS Computational Biology 15:e1006650. https://doi.org/10.1371/journal.pcbi.1006650

Broennimann O, Fitzpatrick MC, Pearman PB, Petitpierre B, Pellissier L, Yoccoz NG, et al. 2012. Measuring ecological niche overlap from occurrence and spatial environmental data. Global Ecology and Biogeography 21:481–497. https://doi.org/10.1111/j.1466-8238.2011.00698.x

Brown JL, and Carnaval AC. 2019. A tale of two niches: methods, concepts, and evolution. Frontiers of Biogeography 11:e44158. https://doi.org/10.21425/F5FBG44158

Burgin CJ, Colella JP, Kahn PL, and Upham NS. 2018. How many species of mammals are there? Journal of Mammalogy 99:1–14. https://doi.org/10.1093/jmammal/gyx147

Burnham KP, and Anderson DR. 2002. Model selection and multimodel inference: a practical information-theoretic approach. New York (USA): Springer.

Calixto-Pérez E, Alarcón-Guerrero J, Ramos-Fernández G, Dias PAD, Rangel-Negrín A, Améndola-Pimenta M, et al. 2018. Integrating expert knowledge and ecological niche models to estimate Mexican primates’ distribution. Primates 59:451–467. https://doi.org/10.1007/s10329-018-0673-8

Carstens BC, Pelletier TA, Reid NM, and Satler JD. 2013. How to fail at species delimitation. Molecular Ecology 22:4369–4383. https://doi.org/10.1111/mec.12413

Castellanos-Morales G, Gámez N, Castillo-Gámez RA, and Eguiarte LE. 2016. Peripatric speciation of an endemic species driven by Pleistocene climate change: the case of the Mexican prairie dog (Cynomys mexicanus). Molecular Phylogenetics and Evolution 94:171–181. https://doi.org/10.1016/j.ympev.2015.08.027

Castellanos-Morales G, Gasca-Pineda J, Ceballos G, and Ortega J. 2014. Genetic variation in a peripheral and declining population of black-tailed prairie dogs (Cynomys ludovicianus) from Mexico. Journal of Mammalogy 95:467–479. https://doi.org/10.1644/12-MAMM-A-099

Cothran EG, Zimmerman EG, and Nadler CF. 1977. Genic differentiation and evolution in the ground squirrel subgenus Ictidomys (genus Spermophilus). Journal of Mammalogy 58:610–622. https://doi.org/10.2307/1380009

D’Elía G, Fabre PH, and Lessa EP. 2019. Rodent systematics in an age of discovery: recent advances and prospects. Journal of Mammalogy 100:852–871. https://doi.org/10.1093/jmammal/gyy179

Dallas TA, and Ten Caten C. 2025. Linking geographic distribution and niche through estimation of niche density. Journal of Animal Ecology 94:1221–1230. https://doi.org/10.1111/1365-2656.70052

Darriba D, Taboada GL, Doallo R, and Posada D. 2012. jModelTest 2: more models, new heuristics and parallel computing. Nature Methods 9:772–772. https://doi.org/10.1038/nmeth.2109

De Queiroz K. 2007. Species concepts and species delimitation. Systematic Biology 56:879–886. https://doi.org/10.1080/10635150701701083

Di Cola V, Broennimann O, Petitpierre B, Breiner FT, d’Amen M, Randin C, et al. 2017. ecospat: an R package to support spatial analyses and modeling of species niches and distributions. Ecography 40:774–787. https://doi.org/10.1111/ecog.02671

Dinerstein E, Olson D, Atchley J, Loucks C, Contreras-Balderas S, Abell R, et al. 2000. Ecoregion-based conservation in the Chihuahuan Desert. Washington, DC (USA): World Wildlife Fund.

Edwards CW, Fulhorst CF, and Bradley RD. 2001. Molecular phylogenetics of the Neotoma albigula species group: further evidence of a paraphyletic assemblage. Journal of Mammalogy 82:267–279. https://doi.org/10.1644/1545-1542(2001)082<0267:MPOTNA>2.0.CO;2

Ezard T, Fujisawa T, and Barraclough TG. 2009. Splits: species’ limits by threshold statistics. R package versión 1:r29.

Fabre PH, Hautier L, Dimitrov D, and Douzery EJ. 2012. A glimpse on the pattern of rodent diversification: a phylogenetic approach. BMC Evolutionary Biology 12:88. https://doi.org/10.1186/1471-2148-12-88

Feijó A, Vilela JF, Cheng J, Schetino MAA, Coimbra RT, Bonvicino CR, et al. 2019. Phylogeny and molecular species delimitation of long-nosed armadillos (Dasypus: Cingulata) supports morphology-based taxonomy. Zoological Journal of the Linnean Society 186:813–825. https://doi.org/10.1093/zoolinnean/zly091

Fernández JA, Cervantes FA, and Hafner MS. 2012. Molecular systematics and biogeography of the Mexican endemic kangaroo rat, Dipodomys phillipsii (Rodentia: Heteromyidae). Journal of Mammalogy 93:560–571. https://doi.org/10.1644/11-MAMM-A-224.1

Fernández JA. 2012. Phylogenetics and biogeography of the microendemic rodent Xerospermophilus perotensis (Perote ground squirrel) in the Oriental Basin of Mexico. Journal of Mammalogy 93:1431–1439. https://doi.org/10.1644/11-MAMM-A-409.1

Ferrari L, López-Martínez M, Aguirre-Díaz G, and Carrasco-Núñez G. 1999. Space-time patterns of Cenozoic arc volcanism in central Mexico: from the Sierra Madre Occidental to the Mexican Volcanic Belt. Geology 27:303–306. https://doi.org/10.1130/0091-7613(1999)027<0303:STPOCA>2.3.CO;2

Ferrari L. 2000. Avances en el conocimiento de la Faja Volcánica Transmexicana durante la última década. Boletín de la Sociedad Geológica Mexicana 53:84–92. http://dx.doi.org/10.18268/BSGM2000v53n1a5

Fišer C, Robinson CT, and Malard F. 2018. Cryptic species as a window into the paradigm shift of the species concept. Molecular Ecology 27:613–635. https://doi.org/10.1111/mec.14486

Flores-Manzanero A, and Vázquez-Domínguez E. ٢٠١٩. Landscape genetics of mammals in American ecosystems. Therya 10:381–393.

Fujisawa T, and Barraclough TG. 2013. Delimiting species using single-locus data and the Generalized Mixed Yule Coalescent approach: a revised method and evaluation on simulated data sets. Systematic Biology 62:707–724. https://doi.org/10.1093/sysbio/syt033

Gamez N, Nihei SS, Scheinvar E, and Morrone JJ. 2017. A temporally dynamic approach for cladistic biogeography and the processes underlying the biogeographic patterns of North American deserts. Journal of Zoological Systematics and Evolutionary Research 55:11–18. https://doi.org/10.1111/jzs.12142

García-Arévalo A, and Nocedal J. 2008. Vegetación de áreas naturales protegidas. In: González-Elizondo SM, González-Elizondo M, and Márquez-Linares MA, editores. Vegetación y ecorregiones de Durango. Durango (MEX): Instituto Politécnico Nacional; p. ١٢١–١٢٨.

GBIF: The Global Biodiversity Information Facility. ٢٠٢٤. What is GBIF? [Consultado el ١٥ de mayo de ٢٠٢٤]. https://www.gbif.org/

Ge DY, Liu X, Lv XF, Zhang ZQ, Xia L, and Yang QS. 2014. Historical biogeography and body form evolution of ground squirrels (Sciuridae: Xerinae). Evolutionary Biology 41:99–114. https://doi.org/10.1007/s11692-013-9250-7

Gearty W. 2025. deeptime: an R package that facilitates highly customizable and reproducible visualizations of data over geological time intervals. Big Earth Data 1–17. https://doi.org/10.1080/20964471.2025.2537516

Gittenberger E. 1991. What about non-adaptive radiation? Biological Journal of the Linnean Society 43:263–272. https://doi.org/10.1111/j.1095-8312.1991.tb00598.x

Goldstein PZ, and DeSalle R. 2011. Integrating DNA barcode data and taxonomic practice: determination, discovery, and description. Bioessays 33:135–147. https://doi.org/10.1002/bies.201000036

González-Ruíz N, Álvarez-Castañeda ST, and Álvarez T. 2005. Distribution, taxonomy, and conservation status of the Perote mouse Peromyscus bullatus (Rodentia: Muridae) in México. Biodiversity and Conservation 14:3423–3436. https://doi.org/10.1007/s10531-004-0548-8

Gould SJ. 2002. The structure of evolutionary theory. Cambridge, MA (USA) and London (UK): Harvard University Press.

Guevara-Chumacero LM, López-Wilchis R, Pedroche FF, and Barriga-Sosa IDLA. 2006. Análisis estructural del gen mitocondrial citocromo y de la región control de Cynomys mexicanus y Spermophilus spilosoma (Rodentia: Sciuridae). Acta Zoológica Mexicana 22:123–125.

Hafner DJ, and Riddle BR. 2005. Mammalian phylogeography and evolutionary history of northern Mexico’s deserts. In: Cartron JL, Ceballos G, and Felger RS, editores. Biodiversity, ecosystems, and conservation in Northern Mexico. Oxford (UK): Oxford University Press; p. 225–245.

Hafner DJ, and Riddle BR. 2011. Boundaries and barriers of North American warm deserts. In: Upchurch P, McGowan J, and Slater CSC, editores. Palaeogeography and palaeobiogeography: Biodiversity in space and time. Florida (USA): CRC Press; p. 75–114.

Hafner DJ, Hafner MS, Hasty GL, Spradling TA, and Demastes JW. 2008. Evolutionary relationships of pocket gophers (Cratogeomys castanops species group) of the Mexican Altiplano. Journal of Mammalogy 89:190–208. https://doi.org/10.1644/07-MAMM-A-001.1

Hansen TF, and Houle D. 2004. Evolvability, stabilizing selection, and the problem of stasis. In: Pigliucci M, and Preston K, editores. Phenotypic integration: Studying the ecology and evolution of complex phenotypes. Oxford (UK): Oxford University Press; p. 130–152.

Hernández-Canchola G, León-Paniagua L, and Esselstyn JA. 2021. Mitochondrial DNA indicates paraphyletic relationships of disjunct populations in the Neotoma mexicana species group. Therya 12:411–421. https://doi.org/10.12933/therya-21-1082

Harrison RG, Bogdanowicz SM, Hoffmann RS, Yensen E, and Sherman PW. 2003. Phylogeny and evolutionary history of the ground squirrels (Rodentia: Marmotinae). Journal of Mammalian Evolution 10:249–276. https://doi.org/10.1023/B:JOMM.0000015105.96065.f0

Helgen KM, Cole FR, Helgen LE, and Wilson DE. 2009. Generic revision in the Holarctic ground squirrel genus Spermophilus. Journal of Mammalogy 90:270–305. https://doi.org/10.1644/07-MAMM-A-309.1

Hey J. 2006. On the failure of modern species concepts. Trends in Ecology and Evolution 21:447–450. https://doi.org/10.1016/j.tree.2006.05.011

Ho SY. 2007. Calibrating molecular estimates of substitution rates and divergence times in birds. Journal of Avian Biology 38:409–414. https://doi.org/10.1111/j.0908-8857.2007.04168.x

Hoang DT, Chernomor O, Von Haeseler A, Minh BQ, and Vinh LS. 2018. UFBoot2: improving the ultrafast bootstrap approximation. Molecular Biology and Evolution 35:518–522. https://doi.org/10.1093/molbev/msx281

Hope AG, Malaney JL, Bell KC, Salazar-Miralles F, Chavez AS, Barber BR, et al. 2016. Revision of widespread red squirrels (genus: Tamiasciurus) highlights the complexity of speciation within North American forests. Molecular Phylogenetics and E´volution 100:70–182. https://doi.org/10.1016/j.ympev.2016.04.014

Jombart T. 2008. Adegenet: an R package for the multivariate analysis of genetic markers. Bioinformatics 24:1403–1405. https://doi.org/10.1093/bioinformatics/btn129

Kass JM, Vilela B, Aiello‐Lammens ME, Muscarella R, Merow C, and Anderson RP. 2018. Wallace: A flexible platform for reproducible modeling of species niches and distributions built for community expansion. Methods in Ecology and Evolution 9:1151–1156. https://doi.org/10.1111/2041-210X.12945

Kelley VC. 1952. Tectonics of the Rio Grande depression of central New Mexico. In: Johnson RB, and Read CB, editores. Rio Grande Country: New Mexico Geological Society Guidebook 3rd. Socorro, NM (USA): New Mexico Geological Society; p. 92–105.

Koprowski JL, Perault DR, Wells MC, and Pasch B. 2016. Family Sciuridae (tree, flying and ground squirrels, chipmunks, marmots and prairie dogs). In: Wilson DE, Lacher TE Jr, and Mittermeier RA, editores. Handbook of the mammals of the world. Vol. 6, Lagomorphs and Rodents I. Barcelona (ESP): Lynx Edicions; p. 648–837.

Kozak KH, and Wiens J. 2006. Does niche conservatism promote speciation? A case study in North American salamanders. Evolution 60:2604–2621. https://doi.org/10.1111/j.0014-3820.2006.tb01893.x

Kulmuni J, Butlin RK, Lucek K, Savolainen V, and Westram AM. 2020. Towards the completion of speciation: the evolution of reproductive isolation beyond the first barriers. Philosophical Transactions of the Royal Society B Biological Sciences 375:20190528. https://doi.org/10.1098/rstb.2019.0528

Kumar S, Stecher G, Li M, Knyaz C, and Tamura K. 2018. MEGA X: molecular evolutionary genetics analysis across computing platforms. Molecular Biology and Evolution 35:1547–1549. https://doi.org/10.1093/molbev/msy096

Lacher T, Timm R, and Álvarez-Castañeda ST. 2016. Xerospermophilus spilosoma. IUCN Red List of Threatened Species V 2018.2. [Consultado el 21 de abril de 2025]. https://www.iucnredlist.org/species/42563/22262899

Leigh JW, Bryant D, and Nakagawa S. 2015. POPART: full-feature software for haplotype network construction. Methods in Ecology and Evolution 6:1110–1116. https://doi.org/10.1111/2041-210X.12410

León-Tapia MÁ, Nava-García E, Cervantes FA, and Arellano E. 2023. Phylogeographic relationships, structure, and genetic diversity of the Mexican endemic volcano Harvest Mouse Reithrodontomys chrysopsis (Rodentia, Cricetidae). Journal of Mammalian Evolution 30(3):735–746. https://doi.org/10.1007/s10914-023-09670-8

Loera I, Ickert‐Bond SM, and Sosa V. 2017. Pleistocene refugia in the Chihuahuan Desert: the phylogeographic and demographic history of the gymnosperm Ephedra compacta. Journal of Biogeography 44:2706–2716. https://doi.org/10.1111/jbi.13064

Lorenzo C, García-Bautista M, Rosas-Ronzón C, Álvarez-Castañeda ST, and Brown DE. 2021. Morphological and genetic variation of black-tailed jackrabbit (Lepus californicus) populations separated by rivers. Therya 12:257–274. https://doi.org/10.12933/therya-21-1104

Losos JB. 2008. Phylogenetic niche conservatism, phylogenetic signal and the relationship between phylogenetic relatedness and ecological similarity among species. Ecology Letters 11:995–1003. https://doi.org/10.1111/j.1461-0248.2008.01229.x

Losos JB. 2011. Convergence, adaptation, and constraint. Evolution 65:1827–1840. https://doi.org/10.1111/j.1558-5646.2011.01289.x

Luna-Aranguré C, and Vázquez-Domínguez E. ٢٠٢٠. Analysis of the application of ecological niche modeling in phylogeographic studies: contributions, challenges, and future. Therya 11:47–55. https://doi.org/10.12933/therya-20-844

Mammal Diversity Database. 2025. Mammal Diversity Database V 1.13. [Consultado el 2 de noviembre de 2025]. https://www.mammaldiversity.org/

Martínez‐Borrego D, Arellano E, González‐Cózatl FX, Ospina‐Garcés SM, and Rogers DS. 2023. Species delimitation and integrative taxonomy of the Reithrodontomys mexicanus (Rodentia: Cricetidae) cryptic complex. Ecology and Evolution 13:e10355. https://doi.org/10.1002/ece3.10355

McLean BS, Rickart EA, Cook JA, Guralnick RP, Burgin CJ, and Lohr K. 2025. Integrative species delimitation reveals an Idaho-endemic ground squirrel, Urocitellus idahoensis (Merriam 1913). Journal of Mammalogy 106:405–429. https://doi.org/10.1093/jmammal/gyae135

Menéndez I, Gomez Cano AR, Cantalapiedra JL, Peláez‐Campomanes P, Álvarez‐Sierra MÁ, and Hernandez Fernandez M. 2021. A multi‐layered approach to the diversification of squirrels. Mammal Review 51:66–81. https://doi.org/10.1111/mam.12215

Morgan PA, and Golombek MP. 1984. Factors controlling the phases and styles of extension in the northern Rio Grande rift. In: Baldridge WS, Dickerson PW, Riecker RE, and Zidek J, editores. Rio Grande Rift: Northern New Mexico. Socorro, NM (USA): New Mexico Geological Society; p. 13–19.

Morin X, and Lechowicz MJ. 2013. Niche breadth and range area in North American trees. Ecography 36:300-312. https://doi.org/10.1111/j.1600-0587.2012.07340.x

Montero-Bagatella SH, and González-Romero A. 2014. Análisis demográfico de una especie en estado crítico de conservación: la ardilla de Perote. Therya 5:737–749. https://doi.org/10.12933/therya-14-249

Moulatlet GM, Merow C, Maitner B, Boyle B, Feng X, Frazier AE, et al. 2025. General laws of biodiversity: Climatic niches predict plant range size and ecological dominance globally. Proceedings of the National Academy of Sciences 122:e2517585122. https://doi.org/10.1073/pnas.2517585122

Muscarella R, Galante PJ, Soley‐Guardia M, Boria RA, Kass JM, Uriarte M, et al. 2014. ENM eval: An R package for conducting spatially independent evaluations and estimating optimal model complexity for Maxent ecological niche models. Methods in Ecology and Evolution 5:1198–1205. https://doi.org/10.1111/2041-210X.12261

Neiswenter SA, and Riddle BR. 2010. Diversification of the Perognathus flavus species group in emerging arid grasslands of western North America. Journal of Mammalogy 91:348–362. https://doi.org/10.1644/09-MAMM-A-102.1

Neiswenter SA, Hafner DJ, Light JE, Cepeda GD, Kinzer KC, Alexander LF, et al. 2019. Phylogeography and taxonomic revision of Nelson’s pocket mouse (Chaetodipus nelsoni). Journal of Mammalogy 100:1847–1864. https://doi.org/10.1093/jmammal/gyz130

Nguyen LT, Schmidt HA, Von Haeseler A, and Minh BQ. 2015. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Molecular Biology and Evolution 32:268–274. https://doi.org/10.1093/molbev/msu300

Nosil P. 2012. Ecological speciation. Oxford (UK): Oxford University Press.

Olson DM, Dinerstein E, Wikramanayake ED, Burgess ND, Powell GV, Underwood EC, et al. 2001. Terrestrial Ecoregions of the World: A New Map of Life on Earth: A new global map of terrestrial ecoregions provides an innovative tool for conserving biodiversity. BioScience 51:933–938. https://doi.org/10.1641/0006-3568(2001)051[0933:TEOTWA]2.0.CO;2

Pahad G, Montgelard C, and Jansen van Vuuren B. 2019. Phylogeography and niche modelling: Reciprocal enlightenment. Mammalia 84:10–25. https://doi.org/10.1515/mammalia-2018-0191

Paradis E, Claude J, and Strimmer K. 2004. APE: analyses of phylogenetics and evolution in R language. Bioinformatics 20:289–290. https://doi.org/10.1093/bioinformatics/btg412

Pečnerová P, and Martínková N. 2012. Evolutionary history of tree squirrels (Rodentia, Sciurini) based on multilocus phylogeny reconstruction. Zoologica Scripta 41:211–219. https://doi.org/10.1111/j.1463-6409.2011.00528.x

Peixoto FP, Villalobos F, and Cianciaruso MV. 2017. Phylogenetic conservatism of climatic niche in bats. Global Ecology and Biogeography 6:1055–1065. https://doi.org/10.1111/geb.12618

Peng S, Ellison AM, and Davis CC. 2025. Incorporating responses of traits to changing climates into species distribution models: a path forward. New Phytologist 248(1):38–49. https://doi.org/10.1111/nph.70402

Petersen MK. 1976. The Rio Nazas as a factor in mammalian distribution in Durango, Mexico. The Southwestern Naturalist 20:495–502. https://doi.org/10.2307/3669866

Peterson AT, Egbert SL, Sánchez-Cordero V, and Price KP. 2000. Geographic analysis of conservation priority: endemic birds and mammals in Veracruz, Mexico. Biological Conservation 93:85–94. https://doi.org/10.1016/S0006-3207(99)00074-9

Phillips SJ, Anderson RP, Dudík M, Schapire RE, and Blair ME. 2017. Opening the black box: An open‐source release of Maxent. Ecography 40:887–893. https://doi.org/10.1111/ecog.03049

Puillandre N, Lambert A, Brouillet S, and Achaz GJME. 2012. ABGD, Automatic Barcode Gap Discovery for primary species delimitation. Molecular Ecology 21:1864–1877. https://doi.org/10.1111/j.1365-294X.2011.05239.x

R Core Team. 2020. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, AT. https://www.R-project.org/

Rambaut A. 2010. FigTree v1. 3.1. Institute of Evolutionary Biology, University of Edinburgh. Edinburgh (UK). https://tree.bio.ed.ac.uk/software/figtree/

Rambaut A, Drummond AJ, Xie D, Baele G, and Suchard MA. 2018. Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Systematic Biology 67:901–904. https://doi.org/10.1093/sysbio/syy032

Raymo ME, and Ruddiman WF. 1992. Tectonic forcing of late Cenozoic climate. Nature 359:117–122. https://doi.org/10.1038/359117a0

Repasch M, Karlstrom K, Heizler M, and Pecha M. 2017. Birth and evolution of the Rio Grande fluvial system in the past 8 Ma: Progressive downward integration and the influence of tectonics, volcanism, and climate. Earth-Science Reviews 168:113–164. https://doi.org/10.1016/j.earscirev.2017.03.003

Riddle BR, Hafner DJ, and Alexander LF. 2000. Phylogeography and systematics of the Peromyscus eremicus species group and the historical biogeography of North American warm regional deserts. Molecular Phylogenetics and Evolution 17:145–160. https://doi.org/10.1006/mpev.2000.0841

Riddle BR. 1995. Molecular biogeography in the pocket mice (Perognathus and Chaetodipus) and grasshopper mice (Onychomys): the late Cenozoic development of a North American aridlands rodent guild. Journal of Mammalogy 76:283–301. https://doi.org/10.2307/1382342

Rocha RG, Leite YLR, Costa LP, and Rojas D. 2016. Independent reversals to terrestriality in squirrels (Rodentia: Sciuridae) support ecologically mediated modes of adaptation. Journal of Evolutionary Biology 29:2471–2479. https://doi.org/10.1111/jeb.12975

Romero-Almaraz MDL, León-Paniagua L, Ramírez-Pulido J, and Vargas C. 2007. Mamíferos pequeños: manual de técnicas de captura, preparación, preservación y estudio. Ciudad de México (MEX): Las Prensas de Ciencias, Facultad de Ciencias, Universidad Nacional Autónoma de México.

Ronquist F, Teslenko M, Van Der Mark P, Ayres DL, Darling A, Höhna S, et al. 2012. MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Systematic Biology 61:539–542. https://doi.org/10.1093/sysbio/sys029

Rozas J, Ferrer-Mata A, Sánchez-DelBarrio JC, Guirao-Rico S, Librado P, Ramos-Onsins SE, et al. 2017. DnaSP 6: DNA sequence polymorphism analysis of large data sets. Molecular Biology and Evolution 34:3299–3302. https://doi.org/10.1093/molbev/msx248

Rundell RJ, and Price TD. 2009. Adaptive radiation, nonadaptive radiation, ecological speciation and nonecological speciation. Trends in Ecology and Evolution 24:394–399. https://doi.org/10.1016/j.tree.2009.02.007

Sánchez-Cordero V, Illoldi-Rangel P, Linaje M, Sarkar S, and Peterson AT. 2005. Deforestation and extant distributions of Mexican endemic mammals. Biological Conservation 126:465–473. https://doi.org/10.1016/j.biocon.2005.06.022

Scheinvar E, Gámez N, Moreno-Letelier A, Aguirre-Planter E, and Eguiarte LE. 2020. Phylogeography of the Chihuahuan Desert: diversification and evolution over the Pleistocene. In: Ramírez-Bautista A, and Hafner DJ, editores. Plant diversity and ecology in the Chihuahuan Desert: Emphasis on the Cuatro Ciénegas Basin. Cham (CHE): Springer International Publishing; p. 19–44.

Schoener TW. 1970. Nonsynchronous spatial overlap of lizards in patchy habitats. Ecology 51:408–418. https://doi.org/10.2307/1935376

Segurado P, Araujo MB, and Kunin WE. 2006. Consequences of spatial autocorrelation for niche‐based models. Journal of Applied Ecology 43:433–444. https://doi.org/10.1111/j.1365-2664.2006.01162.x

Sikes RS, Animal Care and Use Committee of the American Society of Mammalogists. 2016. 2016 Guidelines of the American Society of Mammalogists for the use of wild mammals in research and education. Journal of Mammalogy 97:663–688. https://doi.org/10.1093/jmammal/gyw078

Simonov E, Lopatina NV, Titov SV, Ivanova AD, Brandler OV, Surin VL, et al. 2024. Traditional multilocus phylogeny fails to fully resolve Palearctic ground squirrels (Spermophilus) relationships but reveals a new species endemic to West Siberia. Molecular Phylogenetics and Evolution 195:108057. https://doi.org/10.1016/j.ympev.2024.108057

Simpson GG. 1961. Principles of animal taxonomy. New York (EEUU): Columbia University Press.

Slatyer RA, Hirst M, and Sexton JP. 2013. Niche breadth predicts geographical range size: a general ecological pattern. Ecology Letters 16:1104–1114. https://doi.org/10.1111/ele.12140

Sullivan RM. 1985. Phyletic, biogeographic, and ecologic relationships among montane populations of least chipmunks (Eutamias minimus) in the southwest. Systematic Zoology 34:419–448. https://doi.org/10.1093/sysbio/34.4.419

Triant DA, and DeWoody JA. 2006. Accelerated molecular evolution in Microtus (Rodentia) as assessed via complete mitochondrial genome sequences. Genetica 128:95–108. https://doi.org/10.1007/s10709-005-5538-6

Van Devender TR, and Burgess TL. 1985. Late Pleistocene woodlands in the Bolson de Mapimí: a refugium for the Chihuahuan Desert Biota? Quaternary Research 24:346–353. https://doi.org/10.1016/0033-5894(85)90056-0

Van Tuinen M, O’Keefe K, Ramakrishnan U, and Hadly EA. 2008. Fire and ice: genetic structure of the Uinta ground squirrel (Spermophilus armatus) across the Yellowstone hotspot. Molecular Ecology 17:1776–1788. https://doi.org/10.1111/j.1365-294X.2008.03671.x

Warren DL, Glor RE, and Turelli M. 2008. Environmental niche equivalency versus conservatism: quantitative approaches to niche evolution. Evolution 62:2868–2883. https://doi.org/10.1111/j.1558-5646.2008.00482.x

Waterman JM, Gossmann TI, Brandler O, and Koprowski JL. 2021. Ecological, behavioral and genomic consequences in the rodent family Sciuridae: why are squirrels so diverse? Frontiers in Ecology and Evolution 9:765558. https://doi.org/10.3389/fevo.2021.765558

Wickham H. 2016. ggplot2: elegant graphics for data analysis. Cham (CHE): Springer International Publishing.

Wiens JJ. 2004. What is speciation and how should we study it? The American Naturalist 163:914–923. https://doi.org/10.1086/386552

Wiens JJ, and Graham CH. 2005. Niche conservatism: integrating evolution, ecology, and conservation biology. Annual Review of Ecology, Evolution and Systematics 36:519–539. https://doi.org/10.1146/annurev.ecolsys.36.102803.095431

Wiens JJ, Ackerly DD, Allen AP, Anacker BL, Buckley LB, Cornell HV, et al. 2010. Niche conservatism as an emerging principle in ecology and conservation biology. Ecology Letters 13:1310–1324. https://doi.org/10.1111/j.1461-0248.2010.01515.x

Wiley EO. 1978. The evolutionary species concept reconsidered. Systematic Zoology 27:17–26. https://doi.org/10.2307/2412809

WorldClim. 2024. WorldClim: Global Climate Data. [Consultado el 15 de mayo de 2024]. https://worldclim.org/

Zelditch ML, Li J, Tran LA, and Swiderski DL. 2015. Relationships of diversity, disparity, and their evolutionary rates in squirrels (Sciuridae). Evolution 69:1284–1300. https://doi.org/10.1111/evo.12642

Zhang J, Kapli P, Pavlidis P, and Stamatakis A. 2013. A general species delimitation method with applications to phylogenetic placements. Bioinformatics 29:2869–2876. https://doi.org/10.1093/bioinformatics/btt499

Associated editor: Mónica Díaz

Submitted: September 23, 2025; Reviewed: November 26, 2025

Accepted: December 12, 2025; Published on line: January 30, 2026.

THERYA, 2026, Vol. 17(1):81-98

DOI: 10.12933/therya.2026.6236 ISSN 2007-3364

Figure 1. Inferred phylogenetic tree for the Xerospermophilus spilosoma complex based on mitochondrial cytochrome b gene sequences (690 bp) using the BI and ML methods (a). Circles represent post-bootstrap probability values. (b) Results of species delineation analyses suggesting the presence of between 1 and 4 potential species within the X. spilosoma complex. (c) Haplotype networks observed for the different populations analyzed; lines perpendicular to the main branches indicate the estimated number of evolutionary steps.

Figure 2. Heatmap showing genetic distances (p) for the Xerospermophilus spilosoma complex. Interspecific genetic distances >0.10 are marked in pink, and intraspecific distances ≤0.6 are marked in green, according to the scale. The table presents the genetic distances (p) between the clades and the species close to X. spilosoma.

Table 1. Results of the similarity and ecological niche equivalence analyses and their interpretation. The asterisk (*) indicates a significant difference.

Comparison	Similarity (D) P-value	Equivalences (D) P-value	Interpretation
A vs B	(0.07) 0.45	(0.03) 0.16	Dissimilar, equivalent
A vs C	(0.01) 0.57	(0.005) 0.97	Dissimilar, equivalent
A vs D	(0) 1	(0) 1.00	Dissimilar, equivalent
B vs C	(0.10) 0.12	(0.07) 0.94	Dissimilar, equivalent
B vs D	(0.21) 0.03*	(0.14) 0.88	Similar, equivalent
C vs D	(0) 0.34	(0.06) 0.51	Dissimilar, equivalent

Figure 3. Phylogenetic relationships of the Xerospermophilus spilosoma complex calibrated with temporal divergence estimates. Nodes represent lineage-divergence events, and horizontal bars indicate 95% confidence intervals for estimated dates, scaled to the late Cenozoic era.

Figure 4. Ellipses indicate 95 % confidence areas around the means of each group; comparison between clades: A) A and B; B) A and C; C) A and D; D) B and C; E) B and D; and F) C and D.

Figure 5. Comparisons between the environmental datasets of clades: A) A and B; B) A and C; C) A and D; D) B and C; E) B and D; and F) C and D, using a principal component analysis (PCA), showing the distribution of populations in the space defined by the first two principal components (PC1 and PC2). Bars in the similarity test histograms represent null models, and the red line represents the observed D-value.

Figure 6. Potential distribution map of phylogenetic clades A, B, C, and D generated by ecological niche modeling using the MAXENT algorithm. The colored areas indicate the estimated environmental suitability for each clade, and symbols (triangle, square, rhombus, and circle) represent collection locations.

Appendix 1

Appendix 1. Geographical data of the specimens used and National Center for Biotechnology Information (NCBI) access number. In the species column, X = Xerospermophilus; I = Ictydomys; and C = Cymomys.

Species	Country	State	Latitude	Longitude	Collection ID	NCBI #	References
X spilosoma	USA	Kansas	37.872	-100.964		AF157885	Harrison et al. 2003
X spilosoma	USA	Kansas	35.084	-106.74		AF157911	Harrison et al. 2003
X spilosoma	USA	New Mexico	35.084	-106.738		JX047300	Fernández 2012
X spilosoma	USA	Arizona	32.1080556	-109.555278	ACUNHC 566	PX673950	This study
X spilosoma	USA	Texas	30.0487278	-103.551506	ACUNHC 1376	PX673951	This study
X spilosoma	USA	Texas	31.8174644	-105.688989	ACUNHC 2180	PX673952	This study
X spilosoma	Mexico	Durango	26.523	-104.089		AF157845	Harrison et al. 2003
X spilosoma	Mexico	Durango	26.524	-103.929		AF157846	Harrison et al. 2003
X spilosoma	Mexico	Durango	2953720	624503	AGR01	PX673953	This study
X spilosoma	Mexico	Durango	2950540	623514	AGR03	PX673954	This study
X spilosoma	Mexico	Durango	2953955	624710	UNAM EY1194	PX673955	This study
X spilosoma	Mexico	Durango	2953955	624710	UNAM MVA103	PX673956	This study
X spilosoma	Mexico	San Luis Potosí	24.125	-100.925		DQ106853	Chumacero et al. 2006
X spilosoma	Mexico	San Luis Potosí	24.2	-100.901667		DQ106854	Chumacero et al. 2006
X perotensis	Mexico	Puebla	19.49	-97.489		AF157840	Harrison et al. 2003
X perotensis	Mexico	Puebla	19.49	-97.489		AF157948	Harrison et al. 2003
X perotensis	Mexico	Veracruz	19.587	-97.33		JX047301	Fernández 2012
X perotensis	Mexico	Puebla	19.49	-97.489		JX047302	Fernández 2012
X perotensis	Mexico	Veracruz	19.572	-97.383		JX047303	Fernández 2012
X perotensis	Mexico	Veracruz	32.4166667	-97.8166667	AMS01	PX673957	This study
X perotensis	Mexico	Veracruz	2161070	678473	AMS02	PX673961	This study
X perotensis	Mexico	Veracruz	2161070	678473	AMS03	PX673960	This study
X perotensis	Mexico	Veracruz	2161070	678473	AMS04	PX673958	This study
X perotensis	Mexico	Veracruz	2161070	677627	AMS06	PX673959	This study
I mexicanus	Mexico	Edomex				AF157848	Harrison et al. 2003
C ludovicianus	Mexico	Chihuahua				JQ885590	Castellanos-Morales et al. 2014