Influence of current climate, historical climate stability and topography on species richness and endemism in Mesoamerican geophyte plants

Species richness and weighted endemism based on presence data

Occurrence data for the Mexican geophyte taxa was based on Cuéllar-Martínez & Sosa (2016), comprising 2,629 occurrence points for 507 species (see Table S1 for complete species list and number of occurrences per species; see also Fig. 1 for examples). The monocot geophyte checklist was based on Espejo-Serna (2012), and complemented with recently described species. Georeferences were obtained by: (1) specimens from Mexican herbaria such as: ENCB, IBUG, IEB, MEXU, XAL (acronyms based on Thiers, 2017); (2) biodiversity databases such as the Global Biodiversity Information Facility, GBIF (http://www.gbif.org/), and the Mexican Biodiversity Database, REMIB (http://www.conabio.gob.mx/remib/); (3) complemented by monographs and floras. A grid with a pixel size of 0.5 decimal degrees was overlaid on a map of Mexico (796 pixels in total to cover the complete extent of Mexico). Species richness (SR) and weighted endemism (WE) were estimated directly using the occurrence points considering each species as absent or present on each pixel. SR was then calculated as the total number of species present on each pixel. For endemism we used the measure of weighted endemism proposed by Linder (2001a) and Linder (2001b), in which the area of interest is divided in grid squares labeled by the degrees of latitude and longitude and the occurrence records are mapped on the grids. The index for weighted endemism is divided by the grid-diversity; each species in every grid is weighted by the inverse of its distribution range in a way that a species found in a single grid is scored as 1. Grids with a small number of species in total, but proportionally with many range-restricted species, display more effectively than species-rich grids with proportionally range-restricted species, giving a corrected index of endemism (Linder, 2001a; Linder, 2001b). Weighted endemism is the sum of the reciprocal of the total number of cells in which each species is found. Thus, WE emphasizes areas that have a high proportion of species with restricted ranges (Crisp et al., 2001). SR and WE were calculated in ArcMap 10.3 using the toolbox SDMtoolbox v1.1c (Brown, 2014). We also estimated SR and WE using species distribution modeling (SDM) for the species with more than five spatially unique occurrence data (107 species, see Table S2) as explained below. We decided to use this reduced dataset of species because SDM is known to have usually low performance when modeling species with very restricted distributions especially when using datasets with less than five points.

Species richness and weighted endemism based on species distribution models

We also used ecological niche modeling to predict both species richness (SR) and weighted endemism (WE) using the estimated distributions based on niche models. Ecological niche models were estimated for species that had at least five spatially unique occurrence points at a resolution of 10 min. Model performance was evaluated using the area under the curve (AUC) and the true skill statistic (TSS) for all estimated niche models and we only used the models that showed good model performance (see Supplementary S2 for AUC and TSS statistics). The niche models were estimated with MaxEnt 3.3 (Phillips & Dudík, 2008) using 10 replications and cross validation. Default values were used for other parameters. The environmental inputs for generating these niche models were based on WorldClim v. 1.4 (Hijmans et al., 2005) bioclimate variables at 10 min resolution for average conditions of 1960–1990 (available at: http://www.worldclim.org/version1). These variables were first selected based on a Pearson correlation analysis (selected variables with r values ≤ 0.70) for the climate data extracted from all spatially unique points. The variables used for niche modeling were: Bio1 = Annual mean temperature, Bio3 = Isothermality (Mean diurnal range/Temperature annual range) (*100), Bio4 = Temperature seasonality, Bio12 = Annual precipitation, Bio15 = Precipitation seasonality, and Bio16 = Precipitation of coldest quarter. Continuous suitability maps were converted to binary using a maximum training plus specificity threshold for species presence for each replicate (Liu, White & Newell, 2013). The latter were then summarized on a binary (presence/absence) map for each species modeled, assuming the presence of species in those pixels that were predicted in all replicates and using them to estimate SR and WE at a 0.5 decimal degrees resolution in SDMtoolbox v1.1c (Brown, 2014) with ArcMap 10.3. Additionally, we also calculated species richness and endemism independently for the four most diverse families of geophyte taxa in this study (Amaryllidaceae, Asparagaceae, Iridaceae and Orchidaceae).

Current climate, historical climate stability and topography

Different predictor variables that represent variations in current climate, climate stability from the LGM, and topography were calculated for our study area. Current climate data was taken from the WorldClim database (Hijmans et al., 2005) using Mexico as the geographic extent. For the historical variables of climate stability we calculated the anomalies between averages for the bioclimate variables from the Last Glacial Maximum (LGM) and current climate scenarios. The climate model CCSM was utilized for LGM and it is available at WorldClim (http://www.worldclim.org/version1), we chose this model based on previous research analyzing climate for LGM in Mexico as well as tests comparing performance of the main models for this period (Ramírez-Barahona & Eguiarte, 2013; Varela, Lima-Riberio & Terribile, 2015). The anomalies between current and past climate were calculated as the absolute value of the difference between present and past variables (e.g., the absolute value of the difference between BIO1 at present and average BIO1 during the LGM). Thus, these historical variables reflect how similar the present climate of each pixel is relative to the LGM average, with higher values indicating larger anomalies. For topography we used elevation (EL) and the topographic roughness index (TRI). TRI shows the variance on elevation of the group of adjacent pixels at a given area, thus reflecting areas that have elevated altitude heterogeneity.

Simultaneous autoregressive models

The estimates based on climate and topographic data were used as predictor variables using spatial autoregressive models (SAR) implemented in SAM (Rangel, Diniz-Filho & Bini, 2010). We decided to use spatial regression models because they are suitable for data with spatial autocorrelation (i.e., adjacent pixels are more similar among themselves than expected by chance) as expected for environmental and species distribution data (Kissling & Carl, 2008). Prior to running these analyses all predictor and response variables were log transformed: log 10 (x + 1). The spatial regressions were carried out taking all predictor variables simultaneously and independently using SR and WE based on both approaches (i.e., direct occurrence points and distribution modeling) as response variables. We tested alpha values from one to five units and by visualizing Moran’s I correlograms for all distance units residual errors of spatial correlations were removed. The optimal value of alpha was determined to be five for the SR and WE estimated with SDM and an alpha value of two was used for the dataset based on direct occurrence data. Only the predictor variables that were significant in the regression analyses were used as a measure of their association with the response variable. We utilized the estimated standardized coefficients of these variables as a measure of the extent of their association to explain variation in SR and WE and to discuss the differences in their contribution to the alternative predictor variables according to our research hypotheses.

Geographic patterns of species richness and endemism

The estimation of SR and WE using direct occurrence data and SDM yielded similar results (Figs. 2C–2F), but distribution modeling had more areas with the highest richness and endemism. The estimated SR and WE of Mexican monocot geophytes mostly coincided in the same areas, though the extent of the area varied. Richness was highest in the mountainous areas of central Mexico, on the Trans-Mexican Volcanic Belt, with lower values in the southern areas of the Sierra Madre Occidental and Sierra Madre Oriental mountain ranges and in the Balsas Basin and the Tehuacán Valley (Figs. 2C and 2E). The highest degree of endemism was restricted to the central area of the Trans-Mexican Volcanic Belt and to particular areas in the Tehuacán Valley and the Balsas Basin (Figs. 2D and 2F). Although SR and WE had similar geographic patterns, some important differences were observed in different mountainous regions of Mexico such as the Sierra Madre del Sur, Chiapas and Sierra Madre Oriental that had some of the highest levels of endemism but did not coincide with the highest species richness particularly for the direct estimates based on occurrence data (dark circles Fig. 2B).

Species richness and endemism based on current climate variables for Amaryllidaceae, Asparagaceae, Iridaceae, and Orchidaceae are shown in Fig. S1. Species richness for Amaryllidaceae and Iridaceae had similar geographic patterns with the highest richness in the eastern section of the Trans-Mexican Volcanic Belt decreasing to the western section. For Asparagaceae and Orchidaceae, species richness was greatest in the central area of the Trans-Mexican Volcanic Belt, decreasing in the western section towards the southern area of the Sierra Madre Occidental. The weighted endemism analysis revealed differences in the area occupied by groups, with Asparagaceae and Orchidaceae more restricted. The geographic patterns of weighted endemism were very similar to those of species richness (Fig. S1).

Spatial regression analyses

Spatial regression analyses for the estimates based on direct occurrence data showed that the predictor variables explained SR (R² = 0.35) (Table 1). Variables with the highest standardized coefficients in the regression using SR included the three types of abiotic predictors: current climate, historical climate and topography. These variables were annual mean temperature (MAP), topographic roughness index (TRI) and stability of mean annual precipitation (SMAP) (Table 1). Among them, the variable with the highest coefficient was MAP, which was positively correlated with species richness, suggesting that the areas with the highest current annual precipitation (green shading in Fig. 2D) harbor more geophyte species. SMAP and TRI showed a weaker but significant association, with negative and positive associations with SR, respectively. This suggests that areas that have changed more in their annual precipitation regimes from the LGM to the present (i.e., less stable precipitation, orange shading in Fig. 2B) have higher species richness. Spatial regressions using WE as the response variable showed that the predictor variables were able to explain SR (R² = 0.25). Significant predictors with the highest standardized coefficients were historical and topographic variables (SMAP and ALT). Both predictors were positively associated with WE, suggesting that areas with the highest geophyte endemism are characterized by significant changes in the level of precipitation since the LGM and high elevations (represented in Figs. 2B and 2E by red and dark shading, respectively).

Spatial autoregressive analyses based on estimates of SR and WE using distribution modeling for the species with at least ten occurrence data points (Table 1, Fig. 3) yielded results similar to those of the regressions based on direct occurrence data, though with stronger associations (R² = 0.87 and R² = 0.64 for SR and WE, respectively). SR and WE were predicted by the three abiotic factors: current climate, historical climate and topography (Fig. 3).

Estimation of species richness and weighted endemism based on occurrence data and ecological niche modeling

The availability of online georeferenced data for taxonomic groups and high-resolution environmental information has facilitated research on species richness-environment relationships and has also allowed for an understanding of the historical and current drivers of variation in species richness (Araújo & Guisan, 2006; Araújo et al., 2008; Loarie et al., 2009). However, even comprehensive distribution databases are expected to underestimate species richness. Particularly for plants, inventory is incomplete for many groups and many regions or species have very restricted distributions with a few records and this might affect the reliability of environment-spatial prediction of species richness (e.g., Yan, Yang & Tang, 2013; Meyer, Weigelt & Kreft, 2016). Moreover, it has been demonstrated that larger sample sizes are needed to correctly map the distribution ranges of species (Feeley & Silman, 2011). Many Mesoamerican monocot geophytes have limited distributions and although they have been collected throughout their range they do not reach an elevated number of records, and some of them are known exclusively from the type collection in a single locality. Ecological niche modeling has been proposed as an alternative for improving the estimation of species distribution when sampling might be incomplete for some regions or scarce for certain groups of organisms and we followed this approach (Graham & Hijmans, 2006). Here we found that estimates of SR and WE, whether by occurrence data or by ecological niche models, yielded similar results, although for WE we only included species with wider distributions and those that had good model performance compared with the occurrence dataset that included all species georeferences (Van Proosdij et al., 2016). These results suggest that there is a strong relationship between climate and species richness for geophyte taxa, for both the species that have restricted distributions (e.g., Petronymphe decora, Milla valliflora) and those that have wider distributions (e.g., Nothoscordum bivalve, Sprekelia formossisima, see Fig. 1). Thus, factors limiting the extent of the distribution of some species might be related to dispersal ability, biotic interactions or time since speciation.

Patterns of species richness and weighted endemism

According to our results, the richness and endemism of monocot geophytes are highest in mountain chains such as the Trans-Mexican Volcanic Belt and the southern parts of two mountain ranges: the Sierra Madre Oriental and the Sierra Madre Occidental (see Fig. 2). These patterns were identical when the four most diverse families are analyzed (Orchidaceae: 181 species, Asparagaceae: 114 species, Amaryllidaceae: 106 species and Iridaceae: 92 species; see Fig. S1).

Little is known about climate change in the mountains of Mexico, however it has been reported that during the glacial periods of the Pleistocene temperatures decreased and precipitation regimes changed, as did seasonality along the Trans-Mexican Volcanic Belt (Metcalfe, 2006). In the late Pleistocene the western section of this ridge was wetter than it is at present, while the eastern part was drier than at present (Bradbury, 1997; see also Fig. 3B). At the end of the Pleistocene the Trans-Mexican Volcanic Belt was very dry in many areas and wetter conditions were identified at the beginning of the Holocene (Metcalfe et al., 2000). Seasonal climates during these periods likely favored the establishment of these taxa along this range. Seasonal climates have been identified as suitable for monocots with this life form here and in a previous study (Cuéllar-Martínez & Sosa, 2016). Moreover, the same areas of the Trans-Mexican Volcanic Belt are also associated with a high degree of richness and endemism for other plant groups, such as the perennial species of Cosmos (Asteraceae, Vargas-Amado et al., 2013), tribe Tigridieae (Iridaceae), which includes most of the geophytes studied here (Munguía-Lino et al., 2015), and fern species in the order Polypodiales (Sanginés-Franco et al., 2015).

It has been indicated that rise of temperature is one of the most important drivers of shifts of species ranges in a number of plant life forms such like trees, epiphytes and understorey plants and that species are limited in their capacity to track climate warming (Hsu et al., 2012; Dullinger et al., 2015; Sittaro et al., 2017). Furthermore climate velocity has been calculated in some areas such like North America finding that vectors for climate variables display a complex mosaic of patterns that vary in space and time and are dependent on the spatial resolution of input data (Dobrowski et al., 2013). On the other hand, it has been pointed out that plant movements are complex and that are not realistically represented in models actually utilized to predict future vegetation changes and that many plant species will probably lag behind broad-scale patterns of climate change (Corlett & Westcott, 2013). The grid utilized by us was of 50 × 50 km and based on the historical climate variables and on the predictive ecological niche models, species movements could be estimated. However in view of the complexity of plant movement, the migration lag and the complex pattern of climate change in regions we decided not to estimate this movement. Nevertheless more accurate models in the future might assess the movement of geophytes in Mesoamerica.

Abiotic drivers of variation in species richness and endemism

Our results suggest that although species richness is mostly influenced by the current climate, there is a strong historical component that has shaped the patterns of geophyte species diversity and affected patterns of endemism. Our hypothesis was that species richness would be better predicted by the variability in current climate and that endemism would be better predicted by climate stability from the time of the Last Glacial Maximum to the present and by topographic features. However, spatial regressions showed that historical climate variables also played an important role in shaping species richness, particularly on mean annual precipitation stability. Mean annual precipitation stability was negatively associated to SR suggesting that the areas that currently have higher precipitation respect to the LGM prediction have more species. These results coincide with those reported for other plants and regions. For example, palm species richness in Madagascar was strongly influenced by paleo-precipitation since the LGM, suggesting that long-term climate history played an important role in distribution patterns (Rakotoarinivo et al., 2013) For patterns of species diversity based solely on current climate, instability of variables such as water deficit and temperature has been reported to be associated with differences in plant diversity in different biogeographical regions in China (Zhang, Silk & Ma, 2016). On the Qinghai-Tibetan Plateau Yan, Yang & Tang (2013) found that differences in climate were associated with life form and the richness of woody plants was correlated with all climate variables, while that of herbaceous plants was associated with water availability. Moreover, we suggest that life forms, like the geophytes might be less constrained in their response to climate change and persist in areas with changing climate because of their perennating organs.

Monocot geophytes are clearly a functional group, they are not restricted to a certain region, in contrast they are distributed not only in mediterranean areas but also in zones with seasonal climates in low to high elevations and they belong to different lineages, the majority in the most diverse groups, such like Orchidaceae, Amaryllidaceae, Iridaceae and Asparagaceae. Mountain chains in Mexico have a complex array of habitats and pine or oak forests are found in soils with low drainage where some of the geophytes were recorded (Cuéllar-Martínez & Sosa, 2016). Their underground organs that vary in form, size and position permit growth in poor rocky soils or even in gypsum soils (Cuéllar-Martínez & Sosa, 2016). Geophytes represent an interesting group to study adaptations to seasonal and harsh habitats.

We focused in the role that current and historic climate play in patterns of species richness and endemism, however additional abiotic factors such as soils were not addressed and they probably influence these patterns. In Asparagaceae, some geophytes exclusively grow in gypsum soils in the Chihuahuan Desert and thus further research might analyze adaptations of these geophytes to this type of soil. Moreover, biotic interactions like pollination most probably play an important role as well. Many orchids, the most diverse geophyte group in Mesoamerica, are known to have specific pollination systems and they need to be determined to add understanding to the diversity and endemism patterns found in our study.