To highlight the importance of Mexican scorpion diversity, we compiled taxonomic species lists from 19 countries and one territory following the latest published scorpion systematics (Sharma et al., 2015; Santibáñez-López et al., 2018, 2019; shown in Supplementary Material S1).

To create a preliminary map of the scorpion diversity protected by the Mexican NPA system, we divided Mexico into 25-km2 cells, counted the number of species registered at each cell to find species hotspots, and overlapped our map with the current Mexican NPA map. We decided to keep the 25-km2 resolution criterion given the area size of the NPA system. In this way, it was possible to count the number of species protected by each NPA and at what percentage. This same exercise was carried out to create hotspots considering a higher taxonomic category, such as families. This category was favored because of the number of records. Species records were mainly taken from the Portal de Datos Abiertos UNAM (https://datosabiertos.unam.mx), which harbors data from the Colección Nacional de Arácnidos (CNAN). This database was examined by two of us (CESL and GCF) and all dubious records were discarded. In addition, data from the available literature on Mexican scorpion systematics (e.g. González-Santillán and Prendini, 2013; Santibáñez-López et al., 2016; Contreras-Félix and Francke, 2019) was also considered and included. Records from other databases (e.g. GBIF) were revised but not considered because most of them have species delimitation problems.

The ecological niche-modeling tool helps to identify geographically optimal environmental conditions by relating species records with environmental variables. A bioclimatic profile is created and can be looked for in other geographic areas or other times (Peterson et al., 2011). This approach has been very useful for climate change studies by evaluating how much area with optimal conditions will decrease, increase or shift in the forthcoming decades; information that helps suggesting where money and conservation effort should be focused.

Field records used in this study were taken by national and international scorpion experts, coming mainly from the Colección Nacional de Arácnidos, Instituto de Biología, UNAM which is the most well maintained, corroborated and complete data base for Mexican scorpions. In addition, we included geographical data from the original species descriptions when available. We geographically visualized these data before starting the modeling procedure and analyzed each datum to confirm that the species were observed in these distribution ranges. We also decided to use a minimum record distance of 10-km2 between data as a way of mitigating sampling bias in case it existed in specific parts of the country.

From the scorpion species reported in Mexico, we decided to model the most widely disperse species that passed our filters of quality information (curated by us, with more than five records). Most scorpion species in Mexico have very restricted distributions and consequently few unique geographic records with a 10-km2 resolution (Stockwell and Peterson, 2002). However, following the suggestion made by Pearson et al. (2007), we validated all modeled species with a jackkinfe procedure (all of which had < 25 records) (Table 1) to test the predictive capacity of our models with low data numbers. The jackkinfe procedure proposed by Pearson and collaborators involves performing n-1 replicates, leaving out one single training data each time. Each model is then evaluated to determine whether it suitably predicts the excluded locality. We continued our analyses with the 24 species that performed well in the jackknife procedure (Pearson et al., 2007) (a success rate of at least 60% and P < 0.01).

To choose the ecological niche algorithm to be used in this study, we carried out a pilot test with the five scorpion species with highest recorded numbers (at least 15 records): Centruroides balsasensis Ponce-Saavedra & Francke, 2004, Centruroides baergi Hoffmann, 1932, Chihuahuanus bilineatus (Pocock, 1898), Mesomexovis variegatus (Pocock, 1898) and Paravaejovis spinigerus (Wood, 1863); and modeled their distribution with the nine different algorithms (MaxEnt, GAM, GBM, CTA, ANN, SER, FDA, MARS and RF) provided by the BIOMOD platform (Thuiller et al., 2009) in R (R Development Core Team, 2018). With these five species, record numbers were sufficient to carry out cross-validation, splitting our data in training (70%) and validating points (30%). We then validated our models with the most common validation tests used in ecological niche modeling: the true skill statistic (TSS), the area under the receiver operating characteristic curve (ROC) and the kappa statistic (KAPPA) (Allouche et al., 2006). The algorithms that best performed were MaxEnt and GBM. MaxEnt (Phillips et al., 2006) has already proved to have a good performance for modeling species distribution (Warren and Seifert, 2011). This method is a generative algorithm representing a great advantage when using scarce data (Phillips et al., 2006; Pearson et al., 2007). The GBM (Gradient Boosting Machine) is a machine-learning algorithm that, instead of fitting a model in a parameter space, it fits the model in a function space (Friedman, 2001). The GBM algorithm has shown robustness in species distribution modeling by capturing complex process-environment relationships giving very realistic outcomes (Hjort and Marmion, 2009). To increase MaxEnt performance, it has been suggested to evaluate different parameter combinations (features and regularization multiplier) for each species (Supplementary Material S2) through the ENMeval tool (Muscarella et al., 2014). ENMeval helps identifying the proper features prior to model building. It has been proven that using the default parameters, especially with small samples, can result in very different distribution maps (Muscarella et al., 2014). ENMeval also helps to increase performance when modeling takes place across space and time (Muscarella et al., 2014). The ENMeval tool uses as a first criterion to evaluate MaxEnt models with different features and a regularization multiplier, the sample-corrected Akaike information criterion (AICc). This criterion penalizes model complexity, an element that has to be controlled in order to increase extrapolation capability (Moreno-Amat et al., 2015). Among those models with small AICc, we used two more model-selection criteria: the area under the curve (AUC) value (metric to test predictability) and the parsimony principle. We then created 10,000 random pseudo-absences and used all records as training data.

We used the BIOMOD platform to carry out both algorithms chosen simultaneously. BIOMOD is a modeling system that automatically creates different ensemble types. Ensembles have been helpful in diminishing uncertainty coming from algorithms and climate models (Araújo and New, 2007), being desirable in analyses as the one presented in this work. We kept the mean-weight ensemble rule, which is the weighted sum of the model probabilities when the model performance is at least TSS ≥ 0.7 (Thuiller et al., 2009). To create a binary map for the ensemble, we used TSS maximization as a criterion. We decided to create binary maps to facilitate the interpretation of our results. We also took advantage of the BIOMOD platform to automatically calculate the suitable area that will be lost or gained in the future in the presence or absence of full dispersal (new potential distribution area).

To calibrate our models, we used different geographic areas for each evaluated species considering Mexican ecoregions instead of geopolitical limits (the M from the BAM diagram; Soberón and Nakamura, 2009). With this approach, projections are expected to be more ecologically meaningful (Peterson et al., 2011). Within each specific area for each species, we used monthly climate surfaces to create 19 bioclimatic variables that represent annual, seasonal and extreme climatic trends (Nix, 1986). We are aware that there are other environmental factors that influence scorpion species distribution (see e.g. Polis, 1990; Prendini, 2005); however, in this study we focused on bioclimatic variables since we wanted to evaluate the impact of climate change.

To obtain the bioclimatic variables, we used the biovars function of the dismo package (Hijmans and Elith, 2015) in R (R Core Team, 2017). We then carried out a correlation analysis of climatic variables for each species and discarded variables if their correlation was higher than 0.9, which are not expected to add more information to the model and could even increase overfitting and reduce extrapolation performance (Steen et al., 2017) (Supplementary Material S2). We used a regionalized climatology with a 1-km2 resolution under current (1961–2000) and future (short term, 2015−2039; and long term, 2075−2099) climatic conditions (Centro Meteorológico Nacional http://atlasclimatico.unam.mx/AECC_descargas/) under the worst-case representative concentration pathway (RCP) +8.5 W/m2, and under a more optimistic scenario +4.5 W/m2 (IPCC, 2014). The RCPs are greenhouse gas concentration pathways that relate to how much energy will stay on Earth instead of being projected to space, thus increasing temperature. Consequently, the W/m2 unity stands for the energy in watts that stays on Earth per m2.

We used a regionalized weighted atmospheric circulation scenario (REA, Reliability Ensemble Averaging) of general circulation models (IPCC, 2014), which gives more weight to those models that have better performance in projecting climate in different areas of Mexico (Ibarra-Cardeña et al., 2013). General circulation models are one of the greatest sources of uncertainty in climate-change vulnerability analyses (Steen et al., 2017), and using REA is a way to account for this vulnerability. We also used the ExDet tool (Mesgaran et al., 2014) to evaluate how bioclimatic variables range and combinations change under different climatic scenarios in the future (short and long term) in comparison to the current scenario. This tool allowed us to examine the geographic areas where greater changes will take place and extrapolation will be needed (Supplementary Material S3) (Mesgaran et al., 2014).

Once the species distribution modeling under future climatic conditions were created, we overlapped them with the NPA system in Mexico to evaluate the percentage of each species potential distribution area that will be protected in the future.

Given the current knowledge of the scorpion fauna, Mexico has the highest number of species (294), followed by Venezuela (199) (Fig. 1a; Supplementary Material S1). Furthermore, endemism is high in Mexico with 10.7% (245) of all scorpion species and 9% (19) of the scorpion genera (Supplementary Material S1). Despite having a limited sampling of its territory, Mexico's richness is quite dense (with a density of 0.015 spp/km2, Fig. 1b). Even if higher taxonomical categories are subjective, family diversity is a good indicator of a diversity-rich area (Lourenco et al., 2000). Eight families have been recorded in Mexico (Fig. 1a; Supplementary Material S1), whereas only six can be found in China or Vietnam, both with lower species and genera counts.

Fig. 1.

(a) Scorpion species, genus and family richness in 19 countries and one territory (Southern Africa). In the “x” axis the number of families, number of genera per country is indicated by the intensity of coloration of the circle (lighter represents more genera) and the size of the circle indicates the number of species (the bigger the circle the more species). (b) Correlation of species richness and the area in km2 per country or territory (US is only represented by the area of 11 states in which scorpions have been reported). (c) Heatmap of the number of species per family (in phylogenetic order) and their representation in each country. Warm colors indicate fewer species. Phylogeny was modified from Coddington et al. 2004, Santibáñez-López et al., 2019). (d) A subdadult female of Typhlochactas sp. (Typhlochactidae), a representative from the endemic family in Mexico (highlighted in red in panel c; photo courtesy of Jorge Mendoza).

(1.21MB).

Taken together, this information evidences hotspots of scorpion species diversity in the Mexican states of Baja California Sur, Coahuila, Guerrero, Oaxaca, Morelos and Yucatán (Fig. 2B). These hotspots have between five and nine different species. In particular, Baja California Sur has the greatest species richness, reaching up to nine species in 25 km2, followed by Coahuila and Oaxaca. None of these hotspots is under protection by an NPA (Fig. 2B). Other hotspots with fewer species in Coahuila, Oaxaca, Puebla and Morelos do overlap with existing NPAs: Cuatro Ciénegas, Don Martín, Tehuacán-Cuicatlán and Corredor Biológico Chichinautzin, respectively.

Fig. 2.

Scorpion hotspots and their protection status in Mexico. A) Species records inside a natural protected area, B) species hotspots and C) family hotspots. The area in blue indicates the area of the natural protected area system.

(0.32MB).

In terms of family richness, Baja California Sur was also the state with the highest family richness, reaching up to six different families in a 25-km2 cell (Fig. 2C). Other hotspots with fewer families can be found in Puebla, Oaxaca, Guerrero and Veracruz. Of these family hotspots, only one is protected by an NPA, located in the Tehuacán-Cuicatlán NPA, Puebla state (Fig. 2C).

The range change in bioclimatic variables between scenarios was very small for the short term in both RCPs. Only small areas in the north of the country have changes in ranges in comparison with the current scenario, which reduced the level of extrapolation needed (Supplementary Material S3). There is a similar situation in the long-term RCP 4.5; however, there is a more important extrapolation need in around 30% of the territory for the long-term RCP 8.5. Uncertainty in our ecological niche modeling increases when more extrapolation is needed. Consequently, greater precaution needs to be taken with this RCP (Supplementary Material S3 and S4).

Both individual models, MaxEnt and GBM, as well as the weighted ensembles had a model performance of TSS > 0.7 and full AUC > 0.9. If we assume no dispersal of scorpions in the forthcoming decades, all the evaluated species are expected to reduce their potential distribution area (PDA; see column PDA Loss/ND in Table 1). In other words, geographical regions with optimal climatic conditions inside their current distribution area will decrease in both time periods considered (2015−2039 and 2075−2099). However, if species dispersal ability is considered, new PDAs are expected to appear in other geographic regions, resulting in a less negative net change. If we take into account this net change, in the short term for RCP 4.5 and 8.5 W/m2, 15 out of the 24 evaluated species are expected to have negative net changes. This result means that more species are expected to lose rather than win PDA. Similarly, in the long term this tendency can also be seen, with 18 and 17 species having negative net changes under the two different scenarios evaluated, RCP 4.5 and RCP 8.5 W/m2, respectively (Table 1).

Still, with or without dispersal, the three most vulnerable species are Chihuahuanus labrimanus, Hadrurus concolorus Stahnke, 1969 and Chihuahuanus crassimanus (Pocock, 1898). They lose their entire current distribution area, but also will gain very little suitable area in the future. Consequently, even if they had free dispersal abilities, they would have very few climatically suitable areas to go. Other species that lose their entire current PDA under a future evaluated scenario are Vaejovis rossmanni, Diplocentrum zacatecanus Hoffmann, 1931, Megacorums grubbsi Sissom, 1994 and Kochius sonorae (Williams, 1971). Finally, Vaejovis rossmanni Sissom, 1989, Diplocentrum zacatecanus Hoffmann, 1931 and Chihuahuanus bilineatus were species with very negative net changes.

All modeled species had a percentage of their PDA in an NPA under current climatic conditions, although this varied widely among species and went from a minimum of 0.7 % to a maximum of 55.56 % (Table 2, Supplementary Material S4). The average percentage of protected PDA under current conditions was 12.21 % and it did not differ significantly (Mann-Whitney-Wilcoxon test, P > 0.5) from the PDA that will be protected in the future (short and long terms). Nevertheless, 10 species out of 13 did not reach even 10 % of their PDA inside an NPA under current climatic conditions and up to 17 species did not under future climatic conditions (Table 2).

Positive impacts in terms of PDA caused by climate change could also be seen. The top five species with higher positive net changes were the following species: Diplocentrus reddelli, Mesomexovis subcristatus (Pocock, 1898), Thorellius cristamanus (Pocock, 1898), Thorellius intrepidus (Thorell, 1876), C. ornatus and C. baergi (Table 1). Of these, D. reddelli had over 50% of their PDA protected by a NPA under current conditions and is not being protected in the future (Table 2).