Environmental data.
We used occurrence data of the individuals grouped in the delimited species to perform an ecological analysis. Localities from museum databases (not included in the phylogeny) collected within less than 60 km of an individual with a known genetic identity were also included in the ecological analyses (Fig. 1; Appendix 1). To reduce georeferencing errors, each geographic coordinate was rectified against the known distribution for R. mexicanus sensu lato (Hooper 1952; Hall 1981). We removed duplicate point records and thinned the dataset by setting a distance between localities ≤ 3 km using the R package spThing (Aiello-Lammens et al. 2015). This step was required to reduce spatial autocorrelation and avoid overfitting the model derived from the presence of multiple records in the same 1-km2 pixel. A total of 56 localities were employed in the ecological analysis to test niche differences among delimited species. Individuals from El Salvador (clade IIB; see Results) could not be included because they were collected from a single locality.
Six bioclimatic variables at 1-km2 resolution were downloaded from Wordclim 2.0 (Fick and Hijmans, 2017;http://www.worldclim.org ): Bio1 = Annual Mean Temperature, Bio5 = Max Temperature of Warmest Month, Bio6 = Min Temperature of Coldest Month, Bio12 = Annual Precipitation, Bio13 = Precipitation of Wettest Month, and Bio14 = Precipitation of Driest Month. These bioclimatic variables were selected based on previous studies that highlighted their importance in ecological studies of small mammals (e.g.: Santos et al. 2017; Guevara et al. 2018; Stanchak and Santana 2018; Martínez-Borrego et al. 2022b). Due to the arboreal preferences reported for species of subgenusAporodon (Hooper 1952; González-Cozátl and Arellano 2015), we included in our analyses the variable Vegetation Continuous Field (VCF; Hansen et al. 2002;www .landcover.org ) as a measure of the percentage of vegetation cover. This product was obtained from the MODIS sensor using 24 scenarios corresponding to the year 2020, with an original resolution of 250 m. We rescaled the VCF product to the same resolution as the bioclimatic variables using ArGis 10.8 (ESRI 2020).
Ecological niche modeling (ENM) .
We used MaxEnt (Phillip et al. 2006) to build ENMs for each delimited species using the occurrence points and the environmental variables. To improve the models’ fit and predictive ability, different parameter settings were evaluated (Warren and Seifert 2011; Muscarella et al. 2014) using a “checkerboard” method to partition the training and test data. Different values of the regularization multiplier (from 0.5 to 6, in increments of 0.5) and five combinations of feature classes (L = linear, Q = quadratic, H = hinge, P = product, LQHP) were tested. The optimal combination of parameters for each ENM was estimated in the R package ENMeval (Muscarella et al. 2014) and selected according to the lowest delta AICc (Akaike Information Criterion) score. To obtain the ENMs, calibration areas were defined using the minimum convex polygon. Final models were constructed with 50 replicates, and continuous maps were visualized using the cloglog output. The ENMs were binarized using the 10th-percentile cut-off threshold in ArGis 10.8. Each model’s performance was evaluated using the Partial Roc metric implemented in the Niche Toolbox site (http://shiny.conabio.gob.mx:3838/nichetoolb2/), with 1000 Bootstrap iterations and E=0.05.