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.