Introduction

Artibeus aztecus is a medium-sized phyllostomid bat that inhabits the highlands of Mesoamerica. The three allopatric populations of this taxon are recognized as subspecies (^{Davis 1969}): Artibeus aztecus aztecus, found from Sinaloa and Nuevo León to Oaxaca in México; Artibeus aztecus minor, located from Chiapas, México, to Honduras; and Artibeus aztecus major, found from Costa Rica and Panama.

The subspecies A. a. aztecus was typical in evergreen forests at relatively high elevations in the mountains bordering the Mexican Plateau, as low as 1000 m in cloud forest and as high as 2,400 m in the pine-fir forest, it has been recorded in pine-oak forest, coniferous forest, Abies forest, cloud montane forest, agricultural areas (^{López-González and García-Mendoza 2006}; ^{Segura-Trujillo and Navarro-Pérez 2010}; ^{Briones-Salas et al. 2019}; ^{Cerón-Hernández et al. 2022}); In Veracruz (México) it is considered vulnerable because it inhabits forest fragments but can use riparian vegetation as corridors to cross grasslands (Cerón-Hernández et al. 2022). In the case of A. a. minor, it has been reported in coniferous forest, montane cloud forest, grasslands, areas with secondary vegetation, agricultural landscapes and in human settlements (^{Davis 1969}; ^{Kraker-Castañeda et al. 2017}; ^{Lorenzo et al. 2017}; ^{Medina-Van Berkum et al. 2020}). Artibeus a. major is the only subspecies whose distributional pattern was not associated with conifers, but with “cloud forest” atmospheric conditions (Davis 1969); there are records of the subspecies in tropical premontane rainforest and tropical lower montane rainforest (^{Zamora-Mejías and Rodríguez-Herrera 2017}; ^{Pineda-Lizano and Chaverri 2022}).

Artibeus aztecus is a frugivorous bat. Fruit-eating bats in Artibeus are considered important in seed dispersal (^{Saldaña-Vázquez 2019}), which is essential for forest regeneration and maintenance of plant genetic diversity and composition (^{Wang and Smith 2022}), thereby being crucial to forest conservation and management (^{Jordano et al. 2011}). In central México A. aztecus eats wild figs (Ficus sp.), capuli cherries (Prunus serotine), cypress (Cupressus sp), and Mexican hawthorn (Crataegus Mexicana; ^{Solari et al. 2019}).

Previously, ^{Davis (1969}) treated the three populations as subspecies, having observed only subtle differences in color and some cranial, mandibular, forearm, and phalanx measurements. He also assumed no interbreeding among the three populations. Artibeus a. major is the largest of the three subspecies, and A. a. minor is the smallest, while A. a. aztecus is the least dark subspecies. Later, a study that tested the degree to which the potential distribution of one taxon predicted the geographic distribution of its putative sister taxon and vice versa, using the chi-squared statistic to evaluate statistical significance. The study found that the subspecies A. a. aztecus and A. a. minor have similar ecological niches (^{Peterson et al. 1999}). These conclusions were confirmed with the reanalysis of the data using chi-square test statistic and background similairity test using both I and D metrics (^{Warren et al. 2008}).

As in other groups of vertebrates (^{Fitzpatrick and Turelli 2006}; ^{Zink 2012}; ^{Heinicke et al. 2017}), including bats (^{Roberts 2006}; ^{Datzmann et al. 2010}; ^{Monteiro and Nogueira 2011}; ^{Morales-Martínez et al. 2021}), geographic isolation is likely driving the diversification process between the central and northern subspecies of the A. aztecus distribution. Long-term geographic isolation of populations could lead to the accumulation of genetic or phenotypic differences through neutral or selective processes (^{Baker and Bradley 2006}). If distinct ecological conditions are present in each region, they may stimulate the divergence process (^{Turelli et al. 2001}; ^{Kozak and Wiens 2006}).

The study of the environmental requirements of species and the possible differences between them can be a useful tool in evaluating the taxonomic status of populations (^{Buermann et al. 2008}; ^{Lentz et al. 2008}; ^{Tocchio et al. 2015}; ^{Guevara and Sánchez-Cordero 2018}). Ecological niche-based modeling (ENM) is a tool that permits the exploration of geographic and ecological processes by relating species occurrence records with environmental data (^{Kozak and Wiens 2006}; ^{Phillips et al. 2006}; ^{Kozak et al. 2008)}. ENM may help make taxonomic decisions by making niche comparisons between populations or species or by identifying regions that could isolate them (^{Rissler and Apodaca 2007}; ^{Martínez-Gordillo et al. 2010}; ^{Arribas et al. 2013}; ^{Aguilar 2019}; ^{Hending 2021}).

Here we evaluate the similarities -or differences- between the climatic requirements of the three subspecies of A. aztecus, using background similarity tests (as in ^{Warren et al. 2008}; but we used a higher number of specimens for each subspecies) and comparing its potential geographic distributions to better understand the ecological resemblance of the subspecies and clarify the taxonomic status of this bat across Mesoamerica. Based on previous studies, we hypothesize that niche conservatism has caused the isolation of A. aztecus populations and possible morphological divergence.

Materials and methods

Occurrence data. We collected georeferenced occurrence records for the three populations from the Mammal Collection of the Zoology Museum, UNAM (Facultad de Ciencias - Universidad Nacional Autónoma de México, México City, México, MZFC-M), the Mammal Collection of CIDIIR Durango (Centro Interdisciplinario de Investigación para el Desarrollo Integral Regional, Unidad Durango, Instituto Politécnico Nacional, Durango City, México, CRD), and from the databases of VertNet (downloaded on July 27, 2020) and of the Global Biodiversity Information Facility (GBIF, http://www.gbif.org; downloaded on April 30, 2021: https://doi.org/10.15468/dl.e2b69x), using the name “Artibeus aztecus” recorded from 1960 to the present (2020-2021). To reduce sampling bias, we spatially thinned our original data set using the spThin package (^{Aiello-Lammens et al. 2015}) in R 4.0.3. While retaining the greatest number of localities possible, thinning ensured that the distance between all pairs of localities exceeded 10 km (^{Boria et al. 2014}). Records for the final database are shown in Supplementary material 1. We follow ^{Baker et al. (2016}) and ^{Cirranello et al. (2016}) in using Artibeus rather than Dermanura (contra^{Burgin et al. 2018}).

Environmental data. We used 15 bioclimatic variables (Supplementary material 2; ^{Hijmans et al. 2005}, www.worldclim.org) at ~5 km resolution, excluding the four layers that combine precipitation and temperature information into the same layer since they show odd spatial anomalies between neighboring pixels (^{Escobar et al. 2014}), apparently as a consequence of their linked temperature and precipitation variables (^{Campbell et al. 2015}). We extracted the climatic data using ArcMap (ArcGIS Desktop: Release 10.4).

We used a Pearson correlation test to detect and exclude highly correlated environmental variables. The analysis was performed in R with the library ntbox (^{Osorio-Olvera et al. 2020}), which filters the variables that summarize the environmental information of the presences (occurrences) data according to a correlation threshold; this algorithm suggests which variables to use for the modeling part. The threshold selected for this analysis was r < 0.7.

Calibration area. The dispersal capacity of the species, M of the BAM diagram in distribution theory (^{Soberón and Nakamura 2009}), is useful for choosing the calibration area in niche modeling analysis (^{Barve et al. 2011}). Since the dispersal ability of A. aztecus is unknown, we used ArcMap (ArcGIS Desktop: Release 10.4) to generate the calibration area for each subspecies, with a buffer distance of 1° (~111 km) around occurrences, as a similar distance has been observed in movements of A. lituratus, another species of the genus (^{Arnone et al. 2016}).

Ecological niche modelling. We developed niche models for each of the three subspecies of A. aztecus using the maximum entropy method implemented in Maxent version 3.4.4 (^{Phillips et al. 2006}). To select the models with the optimal settings for each subspecies, we built various models with all the possible combinations of linear, quadratic, and product features, with different percentages of training locations (25 % and 50 %) and different regularization multipliers (from 0.0 to 2.0 in 0.5 steps), analyzing 70 models for each subspecies. We used 10,000 randomly selected pixels within each generated calibration area as the background sample. All the models were generated and evaluated with the library kuenm (^{Cobos et al. 2019}) in R.

We selected the final models based on two evaluation metrics. First, we used partial receiver-operating characteristic (ROC) approaches, as to avoid at least some of the failings of classical ROC approaches (^{Peterson et al. 2008}). We used an acceptable omission error threshold of E = 5 and 100 replicate 50% bootstrap resamplings to establish whether the ROC AUC (area under the curve) ratio was above 1.0. Secondly, we used the 5 % training omission rate (OR05), which shows the proportion of test localities with suitability values lower than those excluding the 5 % of training locations with the lowest predicted suitability. Omission rates above the 10% expectation typically indicate model overfitting (^{Muscarella et al. 2014}). The final models were bootstrapped 10 times and we analyzed the data obtained from the average model.

We analyzed and compared the response curves of the three variables with the highest percentage of contribution and permutation importance for each model. The potential distribution of each subspecies were projected to the Mesoamerican region and to generate binary maps, we chose the 10th percentile training presence threshold (^{Peterson et al. 2007}, 2011). We performed these analyses in ArcMap (ArcGIS Desktop: Release 10.4).

Background similarity test. We used background similarity tests to assess niche differentiation between A. aztecus subspecies (^{Warren et al. 2010}). This test determines whether ENMs are more similar than expected by chance, based on the geographical regions where each subspecies reside. This type of analysis is particularly important when allopatric populations are being compared because some differences in niches may inevitably follow from the fact that distinct geographic regions rarely encompass identical distributions of environmental variables (^{Warren et al. 2010}). We developed 100 replicate comparisons of each population’s known occurrences against the background (points drawn from the accessible area) of the other (sample size matching those available for the “background” population). The background similarity tests were performed with the ENMTools package version 1.0.4 (^{Warren et al. 2021}) in R.

We assess similarity in pairwise combinations of subspecies using two similarity measures: Schoener’s D ⁽¹⁹⁶⁸ and Hellinger’s I. These similarity measures are obtained by comparing the estimates of normalized probability calculated for each grid cell of a study area using a Maxent-generated ENM. Both indexes range from 0, when spaces predicted environmental tolerances do not overlap, to 1, when all grid cells are estimated to be equally suitable for both species. Niche similarity is inferred when the observed value falls above the distribution of expected values. In contrast, the difference is inferred when the observed value falls to the left of the distribution (^{Warren et al. 2010}).

Results

We analyzed 151 confirmed A. aztecus occurrences: 104 for A. a. aztecus, 38 for A. a. minor, and 9 for A. a. major (Figure 1). Ten of the original climate variables were highly correlated with other variables and were excluded from analysis. For the final analysis, we used: annual mean temperature (bio01), mean diurnal range (bio02), isothermality (bio03), annual precipitation (bio12), and precipitation of the driest month (bio14). Final models with the optimal settings for each subspecies were as follow: A. a. aztecus: linear, quadratic, and product features, and regularization multiplier of 1 (Mean AUC ratio: 1.195, OR05: 0.096); A. a. minor: linear and quadratic features, and regularization multiplier of 0.5 (Mean AUC ratio: 1.426, OR05: 0.0); and A. a. major: linear, quadratic and product features and regularization multiplier of 1.5 (Mean AUC ratio: 1.687, OR05: 0.5).

Figure 1 Ocurrence records for the three subspecies of Artibeus aztecus: Artibeus aztecus aztecus (yellow), Artibeus aztecus minor (blue) and Artibeus aztecus major (light red). 

The most important variable for the model of all the subspecies was the annual mean temperature, while the annual precipitation was the only variable that was not placed between the the three most important models for any model. The second and third variable for each model were: mean diurnal range and precipitation of the driest month for A. a. aztecus, precipitation of the driest month and isothermality for A. a. minor, and isothermality and mean diurnal range for A. a. major (Table 1).

Analyzing the response curves for the annual mean temperature, the only important variable in common for the three subspecies, the highest values (> 0.6) of suitability for A. a. aztecus are between 14 °C and 20 °C, while for A. a. minor they are bewteen 14 and 20 °C, and for A. a. major at less at 18 °C (Supplementary material 3). For the mean diurnal range of temperature, in A. a. aztecus the highest suitability is between 8 °C and 14.5 °C, while in A. a. major it is above 6.5 °C (Supplementary material 3a, c). For the isothermality, the highest suitability for A. a. minor was above 70, while for A. a. major it was above 76 (Supplementary material 3b, c). For the precipitation of the driest month, the highest suitability for A. a. aztecus was found at values over 30 mm and for A. a. minor at values between 20 mm and 100 mm (Supplementary material 3a, b).

All potential distribution models showed close correspondence to known distributions of the three populations, showing an association with the highlands of México and Central America (Figure 2). We found relatively wide distributions for the three subspecies, so each model predicted potential distribution areas corresponding with the distribution of the other subspecies. For the three models, the montane regions were separated by less-suitable lowland areas (≤ 500 m), representing potential barriers to the dispersal of each subspecies (e. g., the Isthmus of Tehuantepec and the Nicaraguan Depression).

Pairwise comparisons indicated that A. a. aztecus and A. a. major have the lowest niche overlap (D = 0.246, I = 0.485) and A. a. aztecus and A. a. minor have the highest niche similarity (D = 0.405, I = 0.731). Observed Schoener’s D and Hellinger’s I values were significantly low compared to the null distribution in all cases (Figure 3). Comparisons involving A. a. minor showed D and I values closer to those from the left tail of the null distributions, but significantly different than expected (Figure 3a, c). In sum, background similarity tests indicated that the ecological niche models of the three subspecies were more different than expected by chance (Table 2).

Figure 2 Maxent predicted potential distribution for (a) Artibeus a. aztecus, (b) A. a. minor, and (c) A. a. major.  

Discussion

Potential distributions and geographical barriers. The niche models and potential distribution maps seem to support the findings of the habitat preference of the Aztec fruit-eating bat populations reported previously. Mesoamerican highlands, where the models indicate the potential distribution for each subspecies, include a complex assemblage of montane ecosystems containing high biodiversity and endemism (^{Parra-Olea et al. 2012}; ^{Bryson et al. 2018}; ^{Blair et al. 2019}). Less-suitable areas, such as the Isthmus of Tehuantepec and the Nicaraguan Depression, may act as current geographic barriers to dispersal, limiting contact between the populations, as proposed previously for the subspecies A. a. aztecus and A. a. minor (^{Davis 1969}; ^{Peterson et al. 1999}).

Figure 3 Niche overlap values for Schoener’s D and Hellinger’s I compared to a null distribution: (a) Artibeus a. aztecus (yellow) vs. A. a. minor (blue), (b) A. a. aztecus vs. A. A. major (red), (c) A. a. minor vs. A. a. major.  

Isthmus of Tehuantepec has been proposed as a biogeographic barrier associated with allopatric speciation in a broad range of taxa (^{Sullivan et al. 2000}; ^{León-Paniagua et al. 2007}; ^{Castoe et al. 2009}; ^{Daza et al. 2010}; ^{Rodríguez-Gómez et al. 2013}, 2021) and, climatically, has been considered a barrier for dispersal of oak species, and by separating tropical ecosystems from those with more substantial Nearctic influence (^{Rodríguez-Correa et al. 2015}). The climatic effect of this barrier on the subspecies A. a. aztecus and A. a. minor contrasts with the similar niches found between two haplogroups of the Honduran yellow-shouldered bat Sturnira hondurensis, another Mesoamerican highland bat (^{Hernández-Canchola 2018}).

On the other hand, the Nicaraguan Depression has been considered a major feature determining genetic and biogeographic patterns (^{Gutiérrez-García and Vázquez-Domínguez 2013}). The evolutionary impact of this barrier is reflected in genetic differentiation between sister taxa of vertebrates, including birds (^{Puebla-Olivares et al. 2008}; ^{Arbeláez-Cortés et al. 2010}) and snakes (^{Castoe et al. 2009}). Our findings about the separation between A. a. minor and A. a major are similar to the conclusions of ^{Torres-Morales (2019}), who considered Nicaraguan Depression as a significant barrier that limits the distribution of Sturnira hondurensis, separating it from its sister species S. burtonlimi.

Speciation, and species limits. There is a debate about how conserved the niches between closely related lineages are (^{Wiens and Graham 2005}). Some previous studies have suggested the presence of phylogenetic niche conservatism in phyllostomid bats (Peterson et al. 1999; ^{Stevens 2006}, ²⁰¹¹; ^{Warren et al. 2008}), indicating that closely related species share the same climatic preferences. Alternatively, other authors have not found strong support for niche conservatism in phyllostomid bats (^{Peixoto et al. 2017}), suggesting their niche may have evolved either under strong selection or randomly (^{Diniz-Filho et al. 2010}).

However, former phylogenetic niche conservatism may promote ecological speciation. It can occur in areas with high geographic and ecological variations. In such regions, any geographic distance also results in environmental distance, promoting niche divergence. The combined topographic variation and ecological distance reduce dispersal and gene flow between adjacent populations (^{Gascon et al. 2000}; ^{Gehring et al. 2012}). Lineages may thus adapt to local niches, leading populations to diverge from the ancestral niche (^{Pyron et al. 2015}).

Here, we found signals of ecological niche differentiation among the three subspecies of Aztec fruit-eating bat (Tables 1 and 2, Figures 2 and 3). The three subspecies of A. aztecus present different climatic preferences that may indicate they are evolving independently. Therefore, further studies are necessary to learn about the evolutionary history of A. aztecus and clarify the taxonomic situation of the three subspecies. Certaintly, it is crucial to consider that the outcome and the interpretation of the similarity tests may be sensitive to the definition of the calibration area and environmental background (^{Warren et al. 2010}), still, they may offer some guidelines to explore speciation mechanisms (^{Tocchio et al. 2015}) and thus determine the taxonomic status of the species. In this study, we defined it using the movement data of a congeneric species of A. aztecus, so the results must be carefully interpreted. Further details on the dispersal capacity for each subspecies might improve reference area estimation for niche models.

It is essential to clarify the phylogenetic relationships among the subspecies to better understand their biogeographic history (^{Martínez-Gordillo et al. 2010}). Studies that analyzed the diversification of Artibeus and the subgenus Dermanura, have included a few samples of at least two subspecies, but not A. a. major (^{Owen 1987}; ^{Hoofer et al. 2008}; ^{Redondo et al. 2008}; ^{Solari et al. 2009}; ^{Baker et al. 2016}). Solari et al. (2009) recovered two clades of A. aztecus, represented by samples of A. a. aztecus and A. a. minor, with a genetic divergence of 3.6 % between them, a value that falls in the range necessary for species recognition suggested by Baker and Bradley (2006), so it is crucial to analyze the genetic divergence between the species using a larger number of samples that includes the three subspecies. In addition, morphological analyses that include all subspecies are necessary to assess phenotypic variation and its potential correlation with environmental conditions. A relationship between environmental conditions and morphology has been documented in other Mesoamerican montane species (^{Rodríguez-Gómez et al. 2013}, 2021; ^{Hernández-Canchola 2018}).

In sum, our results offer a first look at the ecological variation of Artibeus aztecus and an additional view on understanding the processes that have shaped the diversification of montane bats in Mesoamerica. Climatic divergence among the three subspecies probably are due to the interaction between former ecological niche conservatism and the emergence of geographic barriers, such as the Isthmus of Tehuantepec and the Nicaraguan Depression that promoted the subsequent ecological differentation.