Material and methods

Occurrence and environmental data. A database with geo-referenced occurrence data of Guadua inermis and Otatea acuminata was built. The information was obtained from field trips, scientific collections herbarium, records from IEB, MEXU and XAL (^{Thiers 2010}) and Global Biodiversity Information Facility (^{GBIF 2015}). Dubious historical records, particularly those with no specimens and those that were duplicates (different collections for the same locality) or offered insufficient geographic information for geo-referencing (i.e., ambiguous localities) were eliminated.

For characterization of the environmental niches we used the 19 bioclimatic variables obtained from the WorldClim Project (^{Hijmans et al. 2005}), with a spatial resolution of 0.0083° (~1 km2) that are the result of interpolating monthly averages from weather stations throughout the world from 1950 to 2000 (Table 1).

Ecological niche modeling and potential species distribution. Prior to the modeling process, we defined polygons to establish the area of accessibility (sensu ^{Soberón 2007}) to run the ENM, since the model does not implicitly take into account the historical barriers of species (^{Barve et al. 2011}). Due to the known distributions of the species, we defined the training area as the accessibility of the species in México and Central America, considering Belize, Guatemala, Honduras and El Salvador. Likewise, we used a digital elevation model (~1 km2/px) for considering the altitudinal range of the species’ distributions.

For modeling performance, we used the algorithm MaxEnt v.3.3.3k (^{Phillips et al. 2006}) which uses the maximum entropy principle to estimate, from the existing values in the climate layers where records occurred, a probability distribution that ranges from 0 to 1 for each pixel, which can be interpreted as an index of habitat suitability for the population that is being modeled (^{Elith et al. 2011}). The algorithm compensates for co-linearity between variables using a method for regularization that deals with feature selection, ranking the contribution of each one through the analysis so there is less need to remove correlated variables (^{Elith et al. 2011}).

For model performance we used for training 80 and 79 % occurrence records for Guadua inermis and Otatea acuminata respectively, and the remaining percentage was used for testing. We used the MaxEnt default values of 500 iterations to limit convergence and 0.00001 convergence limit during modeling, as this is the default level used by background test models. We also used a regularization value of 1. The options for extrapolation and clampling were disabled. Five replicates were set, using a subsample of ~20 % of the records that were run independently to validate each model. To select the best model, we focused on the lowest rate of omission and the highest value of area under the curve (AUC). MaxEnt results indicate the relative suitability of the geographic representation of ecological space in probability values; which were transformed to a binary absence-presence map, using the Minimum Training Presence (MTP) threshold. We chose this threshold considering that all localities used for the model were verified and validated.

The performance of MaxEnt models is traditionally evaluated using the AUC values of the Receiver Operating Characteristic (ROC) curve (^{Phillips et al. 2006}) which allows the evaluation of the coincidence of the climatic suitability generated by the model with the known occurrences, ranging from zero to one, where one indicates perfect discrimination between presences and absences (presence/background in the case of MaxEnt) and 0.5 indicates that the discrimination is no better than the suitability given by a random assumption (^{Fielding & Bell 1997}). However, several problems have been associated with this approach (^{Peterson et al. 2008}, ^{Lobo et al. 2008}), for example, that the two error components (omission and commission) are inappropriately weighted equally. Therefore, we used the partial ROC approach that solves this problem by evaluating only the portion of the ROC curve that covers the spectrum of the prediction, allowing a differential weighting of the two error components (^{Peterson et al. 2008}). Within a value range from 0 to 2, values over 1 suggest a better performance than chance, by analyzing the presences versus the absence against the total area predicted by MaxEnt (^{Barve 2008}). Because of this, we calculated for each model the partial ROC test using the Tool for Partial-ROC v. 1.0. (^{Barve 2008}). We withheld ~20 % of the occurrence records to validate each model. AUCs were limited to the proportional areas over which models actually made predictions and only omission errors < 5 % were considered (^{Peterson et al. 2008}).

Based on the current potential distribution, we projected into future layers under two Representative Concentration Pathways (RCPs) scenarios, which are described in the 5th Assessment Report (^{IPCC 2013}) as externally imposed perturbation in the radiative energy budget of the Earth´s climate system (^{Pachauri et al. 2014}). We used the medium and high emissions projections (RCP 4.5 and RCP 8.5 respectively) in a median future for the year 2050. We used the Centre National de Recherches Météorologiques, France (CNRM) and the Canadian Centre for Climate Modeling and Analysis (CCCma) based on the evaluation of several general circulation models (GCM) by ^{Fernández-Eguiarte et al. (2015)} considering a better characterization of Mexico for the mid future according to the observed data. To reduce uncertainty, we then constructed an ensemble map of both GCM outputs (CCCma, CNRM) for the current and future conditions in order to balance accuracy and robustness by obtaining the coincident areas of geographic potential distribution projected for both models. An ensemble is a collection of climate scenarios where the individual scenarios are different from each other. It highlights the relatively reliability where the different scenarios give similar results than in different directions independently (^{Zhang et al. 2015}).

Results

After a careful revision, we gathered a total of 47 occurrence records for Guadua inermis and 117 unique occurrence records (latitude–longitude) for Otatea acuminata. These occurrence numbers allowed us to develop a satisfactory model performance (high AUC values and significant ratio values of partial ROC analysis; Table 2).

The most important variables identified by the percent of contribution tests (Figure 3) showed that the trend of contribution for each variable is similar in both Canadian and French scenarios. The mean diurnal range (BIO2) and the precipitation in the wettest month (BIO 13) have the highest contribution for G. inermis (Figure 3a), while the temperature and precipitation seasonality (BIO4 and BIO15 respectively) have the highest contribution for O. acuminata in both projections (Figure 3b).

Figure 3 Percent of relative contribution of the 19 climate variables extracted from the analysis. a) Guadua inermis, b) Otatea acuminata. A squares = CCCma; B circles = CNRM. 

To describe bi-dimensionally the ecological niches for both species, we analyzed the relationship between the annual mean temperature and annual precipitation for the occurrence data, and we found that the known distribution of G. inermis occurs in a range of mean annual temperature between 20.3 and 27.2 °C, and annual precipitation between 882 mm and 2,455 mm (Figure 4A). Regarding O. acuminata, the mean annual temperature range goes from 12.5 to 27.8 °C and the annual precipitation from 437 to 1,610 mm (Figure 4B).

Figure 4 Relationship between annual precipitation and annual mean temperature; A) Guadua inermis, B) Otatea acuminata. 

The potential current distribution map revealed that G. inermis (Figure 5A) should be found in the coastal plain of the Gulf of Mexico and the western area of the Yucatan Peninsula under 1,000 m a.s.l, with the exception of the southern part of the distribution in the Chiapas-Guatemala plateau, where it rises to 1,900 m a.s.l. Small areas of the Pacific coastal plains were projected, with attributes similar to those where the distribution of G. inermis is known in the coastal plains of the Gulf of Mexico.

Figure 5 Assemble maps of both GCM outputs (CCCma, CNRM). Green = coincidence of projected distribution with CCCma and CNRM: Blue = projected distribution with CNRM only; Red = projected distribution with CCCma only. A) Guadua inermis, current geographic potential distribution, black dots are geo-referenced localities used to build the models. B) G. inermis 2050 projection, RCP 4.5. C) G. inermis 2050 projection, RCP 8.5. D) Otatea acuminata, current geographic potential distribution, black dots are geo-referenced localities used to build the models. E) O. acuminata 2050 projection, RCP 4.5. F) O. acuminata 2050 projection, RCP 8.5. 

The potential current distribution map for O. acuminata (Figure 5D) suggests that it should be widely present in the Pacific slope, the Balsas depression, the Sierra Madre del Sur, the Sierra Madre Occidental, the Trans-Mexican Volcanic Belt, in some fractions of the coastal plain of the Gulf of México, the Sierra Madre Oriental and in the Central depression in Chiapas. This species has an altitudinal distribution range from 150 to 2,000 m a.s.l. from the records used to calibrate the model.

Because the different model projections do not always correspond to the limits of the predicted potential distribution for each species, the net changes of the ensemble are different to changes in each model independently, namely compiling is spatially explicit the coincident geographic areas between the two models seeking greater certainty in areas predicted (Table 3). We performed an assemble maps in order to find areas of concordance spatially explicit for both GCM, the maps of potential distribution of both species shrink (Figures 5B, C, F, G, Table 3). The changes projected in the assemble map were greater for G. inermis as it loses nearly half of the current distribution (Figure 5C) under RCP 8.5. However, when each of the projections is analyzed an increase of area for the CNRM scenario is shown, in contrast with the Canadian scenario, principally in G. inermis.

Discussion

The analysis of potential effects of climate change on biodiversity via ENM represents a satisfactory approach to understand the tendencies of species distributions in the future (^{Wiens et al. 2009}). In this sense, it is particularly important to pay attention to the impacts of climate change on those species that are endangered, ecologically or geographically restricted, but also those with human economic importance or potential use (^{Stern 2007}).

The studies focused on the future projections of bamboo species in Northern China highlighted the declining trend in distribution (^{Tuanmu et al. 2013}, ^{Li et al. 2015}). This is to be expected for temperate woody bamboos as they are more vulnerable to increases in global temperature, and thus the areas with optimal conditions for them are likely to decrease, whereas those of tropical species could be expected to increase. Furthermore, the output differences from the use of different GCM scenarios generate high uncertainty, and this is why an assemble map is recommended for showing the coincidence of different models and thus, less uncertainty of interpreting results for those areas (^{Zhang et al. 2015}). The species in this study are distributed in a biogeographical transitional zone (Neartic vs Neotropical regions); however, the variables identified with the major contribution of each species suggest different ecological niches for each species. Guadua inermis grows in tropical sub-deciduous forests with less variation in temperature and more annual precipitation, contrary to Otatea acuminata, which grows in tropical dry forests or xerophytic scrub with a large variation in temperature tolerated, especially higher temperatures and less precipitation than G. inermis (Figure 3A, B).

The differences between these contrasting ecological niches are well reflected in the future projections based on assemble maps. For G. inermis the predicted reduction in distribution ensemble is 9.5 and 42.3 % for RCP 4.5 and RCP 8.5 respectively of the total area that was projected for both GCM relative to its current distribution (Figure 5B, C). For O. acuminata the reduction of potential distribution ensemble ranges between 14 to 22 % for RCP 4.5 and RCP 8.5 respectively of the current distribution (Figure 5F, G). Overall for both species, most of this uncertainty will tend to occur on the Pacific slope, the Balsas depression and the Gulf of Mexico coastal plains (Figure 5B-C). However, new areas might be gained mainly in the western portion of the Sierra Madre Occidental and in the northern and southern portions of the Sierra Madre Oriental. These findings do not suggest a threat to these species, and their use and potential management may continue for a long time, independent of climate change. Their management as a provider of culms for wood (for construction, charcoal, furniture, handcrafts) should follow guidelines similar to other commercial species of bamboo, and this management may benefit from anticipating potential areas for cultivation in the following decades. We found climate stability by an ensemble of the maps of current and future projections (both scenarios) for G. inermis on the coast of the Gulf of Mexico and the east of the Sierra Madre del Sur and some parts of the west coast and south of the Yucatan Peninsula (Figure 5D). For O. acuminata the areas of climate stability occur in the Sierra Madre Occidental, the Trans-Mexican Volcanic Belt and the Sierra Madre del Sur and less area in the Sierra Madre Oriental (Figure 5H).

The CC scenarios predicts rising in global temperatures for 2,100 between 1.4 to 3.1 °C under RCP 4.5 and between 2.6 to 4.8 °C under RCP 8.5, and change in the precipitation averages between 10 and 20 % less than today for Mexico (^{Pachauri et al. 2014}). The question is: why would we expect reduction in distribution areas for two species found in fairly high temperature ranges? In the case of G. inermis, a species currently found mostly in lowlands, although temperatures at higher altitudes may increase under CC scenarios, they may not be matching the required increases in precipitation. Therefore, this results in a reduction in the suitable area for large tropical woody bamboos, known to require significant water availability, especially in the season of shoot growth that starts during the rainy season (May or June) (^{Banik 2015}). Otatea acuminata, a species widely distributed in Mexico, seems to be more tolerant of drier climates (^{Ruiz-Sanchez 2015}) and to a little bit lower temperatures than G. inermis. This combination of tolerance climatic variables may make this species more adaptable to climate change and less prone to reduction in geographical distribution compared with G. inermis, especially under the RCP 8.5 scenario. The maintenance of suitable areas for G. inermis and O. acuminata under scenarios of future climate change could keep them playing an important role in the dynamics of the woody vegetation both species inhabit.

Conclusions

As with many clump forming commercial species, both native species are easily propagated, either from vegetative cuttings (rhizomes, culms or branches) or fruits. The wider range of environmental conditions in which Otatea acuminata is found, along with its tolerance to lower temperature and precipitation, should be considered in future incentives to establish commercial plantations and to look for a wider range of traditional and technical uses of this species. In either case, the use of genotypes obtained from sites that more closely match the environmental conditions of the plantation should provide the best results in terms of product quality. It is possible to find some culms of Guadua inermis flowering every year and to get fruits that could be planted for seedlings. On the other hand, although O. acuminata is said to flower every 30 years (^{Ruiz-Sanchez et al. 2011}), it is usually possible to find at least one population flowering every year in different regions in Mexico due to the asynchronic flowering cycles of this species (^{Ruiz-Sanchez et al. 2011}). Both species can be grown from seeds to adults in only about seven years, when they reach the final size to start harvesting culms to be used in construction or to satisfy the growing demand nationwide and worldwide as a building material, as well as handcrafts, baskets and other many potential uses. Promoting their cultivation could help to reduce some climate change effects, for example using both species to restore landslides on hillsides and for CO2 sequestration.