Projections for Mexico’s Tropical Rainforests Considering Ecological Niche and Climate Change

: The tropical rainforest is one of the lushest and most important plant communities in Mexico’s tropical regions, yet its potential distribution has not been studied in current and future climate conditions. The aim of this paper was to propose priority areas for conservation based on ecological niche and species distribution modeling of 22 species with the greatest ecological importance at the climax stage. Geographic records were correlated with bioclimatic temperature and precipitation variables using Maxent and Kuenm software for each species. The best Maxent models were chosen based on statistical signiﬁcance, complexity and predictive power, and current potential distributions were obtained from these models. Future potential distributions were projected with two climate change scenarios: HADGEM2_ES and GFDL_CM3 models and RCP 8.5 W/m 2 by 2075–2099. All potential distributions for each scenario were then assembled for further analysis. We found that 14 tropical rainforest species have the potential for distribution in 97.4% of the landscape currently occupied by climax vegetation (0.6% of the country). Both climate change scenarios showed a 3.5% reduction in their potential distribution and possible displacement to higher elevation regions. Areas are proposed for tropical rainforest conservation where suitable bioclimatic conditions are expected to prevail.


Introduction
Tropical ecosystems are essential because of the ecosystem services it provides (food, medicines), its biological richness and its relevance in supporting life [1]. Tropical rainforest (TR) are plant communities found in warm humid areas and include tropical rainforest and tropical sub-evergreen forest communities. Their dominant species have a perennial physiognomy, that is, they are trees that are more than 30 m high and do not change their leaves seasonally [2]. It is the richest vegetation in floristic composition and biotypes in Mexico's tropical areas, but it currently has a limited and fragmented distribution. In addition, much of it is in early successional stages [3]. In the country, TRs in a mature state (old growth or climax) occupy 12,250 km 2 (0.6% of the national territory) according to the last official report in 2016 [4]. The rest of its successional stages occupy 1% of the country's area; the original distribution area is unknown.
There is greater anthropogenic pressure on TRs than on temperate forests and they are subject to greater landscape fragmentation [5]. During the last century in Mexico, at least 90% was cut down and its use was changed to agriculture, livestock or agroforestry purposes marked by low floristic diversity [3]. Landscapes with substantial differences in plant cover [6] or landscapes fragmentated by deforestation can be found in Veracruz or Chiapas southern states [7]. In addition, climate change has recently been added as a relevant challenge for the ecosystem [8,9].
Mexico's location between two oceans and its complex relief lead to it being particularly impacted by climate change. For the period 2015-2039, annual temperature increases of 1-1.5 • C are estimated for most of the country and up to 2 • C for the northern zone. For precipitation, decreases of 10-20% are projected, as well as variation in its seasonal distribution. El Niño variability on regional scales likely intensifies [10]. Given these circumstances, it is necessary to know how the change in climate will affect the future geographic distribution of TRs in order to promote actions in the short term aimed at their conservation. There is evidence that many terrestrial species have shifted their geographic ranges in response to ongoing climate change [11]. TRs support great associated biological diversity [1], which may be indirectly affected by climate change. To assess the effects of climate change on natural systems such as ecosystems, long-range time horizons are recommended, such as 2075-2099 [12,13]. With the above, possible adverse effects on the level of natural systems are covered, since ignoring long periods could lead to an underestimation of their impact. In addition, a vision at a higher organizational level than the specific one (for example, of plant community or ecosystem) is needed to assess the possible impacts that climate change can have on the potential distribution (PD) of most of their floristic components [14,15]. In this study, the potential distribution area is defined as the land occupied by a species plus invadable distributional area, where the abiotic and biotic conditions are suitable (see Section 2.3 below).
There are several methods for modeling the geographic and potential distribution of species. Some use environmental information as a dependent variable and species records as independent variables. These can be counts, presences and absences, or only presences [16]. The most common method is Ecological Niche Modeling (ENM) and associated Species Distribution Modeling (SDM) [17,18], which have their theoretical bases in the Hutchinson duality [19]. This refers to the environmental conditions where the species lives (ENM) and its corresponding geographic space (SDM) [17,18]. Among the best known techniques and software are BIOCLIM, Generic Algorithm for Rule-Set Prediction (GARP), Generalized Linear Models (GLM), Artificial Neural Network (ANN) and Maximum Entropy (Maxent) [19][20][21][22]. Results obtained with the different techniques may be subtly or dramatically different and may vary due to differences in the performance of different software [21,23]. In the case of Mexico and TR communities, the PD of some species has been previously studied using ENM and SDM with the GARP method [1]. However, higher performance has been obtained with Maxent [22,24], which was selected in this study.
The aims of the work were: (1) to obtain the potential distribution of 22 species that have ecological value and are representative of tropical rainforests under current conditions and future climate change ones, by applying the ENM and SDM techniques through Maxent and Kuenm-they were chosen for their good performance, use and operation [25]; (2) to intersect the areas with the greatest potential distribution in the current and climate change scenarios, in order to identify areas of potential change (increase, decrease, permanence) of the ecosystem; (3) based on the previous objective, to propose areas for conservation of TRs based on climate change and the distribution of their species, also considering current land use and official protection areas. The results contribute to the debate in two ways: by modeling, with cartographic techniques, the potential distribution of an ecosystem based on its species, as other studies have proposed [9,26], and by discussing and providing basic information for decision-makers on environmental legislation and protection. Anticipating the effects of climate change should serve as a tool for proposing areas for conservation based on the richness of potential distributions.

Species Selection
The species were selected for their ecological importance according to the old growth or climax stage of the TR (literature review) and also for their structural value. The ecological importance was based on information extracted from the works of Miranda and Hernández [2], Pennington and Sarukhán [1] and Rzedowski [27]. All of them are Mexican biological and plant species references. Afterwards, two structural valuation indices were applied: importance value (IVI) [28] and forest value (FVI) [29]. Both are common in ecological studies [30]. IVI evaluates dominance, density and frequency, while FVI evaluates diameter at breast height, height and crown cover of the species. Information for the indices was obtained from the National Forest and Soil Inventory (NFSI) for the period 2004-2007 [31,32] and from [4], spatially managed in ArcMap (ESRI) (version 10.5). The results allowed the selection of 22 species (Table 1); the ranking of the ecological and forest value of the species is presented in Supplementary Material S1 (SM1).

Geographic Records and Bioclimatic Variables
The presence data for each species or geographic records (GRs) were obtained from the NFSI [31] and the Global Information Facility [38]. Prior to downloading, GBIF information was filtered by country and geographic coordinates. The information was reviewed to eliminate duplicate data as well as those outside the natural range of the species (urban areas, lakes or oceans). To reduce the sampling bias that often occurs with this type of data (presence only), their spatial autocorrelation was reduced, as recommended by some authors [39], by random filtering or thinning in 5 or 10 km grids. An alternative to using grids is to use the spThin R package [39]. A mesh spacing of 5 km was used when the number of GRs was low (<1500 GR) and 10 km when it was high (>1500 GR). The databases were divided into three spatially unrelated groups with ENMeval, an R package developed by Muscarella et al. [40], to be used in the following modeling phases: training (57 to 69% GR), testing (20 to 32% GR) and final evaluation (10 to 12% GR) of the models. These percentages were varied in order to ensure good representativeness of the group used for the training of the models.
Regarding current bioclimatic variables, the 19 variables developed for BIOCLIM [21] in 1996 and described by O'Donnell and Ignizio [41] were considered. These are biologically significant when combining precipitation and temperature over the course of the year and not only the annual averages [42]. The bioclimatic variables of the current scenario were obtained from WorldClim [43], after which their level of detail was increased according to Gómez et al. [44], who suggest following the delimitation of areas with homogeneous temperature and precipitation called climatic influence areas. Two climate change scenarios were applied: HADGEM2_ES and GFDL_CM3, both for RCP 8.5 scenarios to 2075-2099. The use of these climate change scenarios has been appropriate for climate change studies [9,36] and is recommended since the uncertainty in natural systems is considerable [12]. The General Circulation Models (GCM) were obtained from Fernández et al. [45]. The applied scenarios are considered to represent an increase of up to 4 • C in the average temperature for that period [46], which may be too great for some ecosystems to cope with [11]. The climate change scenarios were also processed using the methodologies of [41,44] to be comparable with current scenarios.

Niche and Species Distribution Modeling
Area M or the accessible area of the BAM diagram by Soberón and Peterson [47] was considered to limit the PD of the species to the space that they could colonize when the interactions between them are omitted [18,47]. Area M refers to the ecological and geographic space where species have the potential to disperse, either because of their historical or relatively recent condition [18,19]. Thus, six M areas were defined for equal groups of species with differences in their distribution. For the above, the country's terrestrial ecoregions were taken as a reference [48] as recommended by [49]. It was then fitted according to the geographic distribution reported in the literature as well as that presented by the GRs for each group of species. The bioclimatic variables were spatially cut to fit the M areas.
A pre-selection of variables was made based on Pearson's correlation coefficient in order to avoid overfitting the models [20] by including variables that are not biologically important [49] and that are redundant in terms of the information they provide. With Maxent, the selection of models free of overfitting is not a problem as it can deal with collinearity between environmental covariates with the regularization choices [49].
Both ENM and SDM were obtained with Maxent (American Museum of Natural History, New York, NY, USA. version 3.4.0) [25] and with the Kuenm R package (The University of Kansas, Lawrence, KS, USA version 1.1.5) [50] in RStudio [26]. The ENM stages consisted of the creation and evaluation of candidate models, and the creation and final evaluation of the best models. The candidate models were created by correlating the pre-selected variables with the GRs and by varying the parameterizations of the Maxent models (regularization parameters and features). The evaluation of the candidate models for the selection of the best model was based on statistical significance with the partial area under the Receiver Operating Characteristic curve (partial ROC) and predictive power (performance) with 5% omission error rates. Complexity was estimated with the Akaike Information Criteria corrected by sample (AICc) according to [51] and implemented in the Kuenm package [50].
The partial ROC criterion was used to verify that the model is not equal to one given by chance by presenting a non-significant p-value [52]. To the extent that the partial ROC tends to 0, there is greater significance in the model. In fact, this statistic is more appropriate than the full area under curve (AUC) for assessing statistical significance [52]; for details of the calculation, see [50]. With the omission rate criterion, the inability of the model (created with training GRs) is obtained to correctly predict the test GRs at a permissible error [53]. The original AIC statistic provides an estimate of the relative information that is lost when varying the models, since it is a measure of compensation between their goodness of fit and their complexity [54]. Thus, complex models are penalized and preference is given to simpler models [55].
The best potential distribution models were created with the training and test GRs as well as the optimal parametrizations using 10 bootstrap replicas with logistic output. The logistic output was interpreted as a level of climatic suitability for the species [9,49]. The best models were evaluated with the final evaluation data with the partial ROC and omission rate criteria using the same package (Kuenm).

Mapping Potential Distribution and Percent Changes
The binary potential distribution maps of each species were reclassified into nonpotential distribution (NPD) and potential distribution (PD) categories. PD was defined with values equal to or greater than 0.3, while NPD refers to values lower than 0.3. Thus, a binary distribution layer was obtained for each species and for each scenario. Using cartographic processes, changes in PD were identified by applying climate change scenarios for each species. The PD in the current scenario was compared to the PD of both future scenarios through the geographic unions of the maps. Additionally, the current PD of each species was compared with the PD obtained in a previous study [1] to contrast the differences shown by the GARP and Maxent methods, as well as by GRs (until 2005 and 2019, respectively).
Plant communities necessarily include ecological communities, since the former provide the necessary conditions for the latter to develop [56,57]. Thus, our idea is that the trends in PD of the main species can indicate the long-term trend of the plant community (tropical rainforest). For this reason, the PDs of all species in each scenario were assembled, grouping them into three potential distribution classes but now with the TR according to the accumulated percentage. Adequate-PD was defined as occurring when a coincidence of 14 to 22 species was found, Possible-PD when there were 8 to 13 species and Non-PD when there were fewer than 7 species. This information was used for the conservation area proposal. This is based on the perspective of the integrated plant community theory, which indicates that where favorable environmental conditions for a plant community are repeated, the species will reappear. This is because the classification is based on similarities [58,59]. Plant communities lack precise borders [2,56], so our hypothesis is based on the idea that the joint presence of species approximates the TR area. Some species require the presence of others in order to coexist, either at the same time or at different times [56]. Regardless of the species, the interactions between species are key in determining the ecological niche of a single species and, therefore, its range [47].

Proposed Conservation Areas
In the proposal for TR conservation areas, priority was given to areas where, when considering climate change, future PDs are expected to prevail or increase. Nationals, state, and voluntarily-conserved protected natural areas [60,61] as well as land use and vegetation [4] were compared to areas where the greatest similarity in prevalence or increase in future PDs is projected by climate change scenarios.

Geographic Records and Bioclimatic Variables
Some authors point out that the number and quality of GRs directly influence the performance of the models, and consequently, the approximation of the species' PD [42,49]. In the present study, an adequate number of GRs was used, since after 14 cleaning and thinning the lowest number of GRs obtained for a species was 85, and Maxent performs well with less than 25 GRs [49,62]. We can see that sampling errors in large databases (such as GBIF) are common due to the heterogeneity of the conditions in which the data are obtained [17], which leads to sampling biases [24,49] and, of course, to duplicate records. Sampling biases also occur because certain areas are not intensively sampled or are not even sampled at all. In other cases they are oversampled [17]. The location of the samples is influenced by accessibility and sometimes by weather conditions, causing the probability of the SDM to favor the oversampled areas [49]. When GR thinning is done, this sampling bias is reduced in areas that have been oversampled but the risk of omitting important data is increased (due to lacking information) [62]. For this reason, the thinning procedure was applied to GRs in a random manner to avoid systematic software errors.
The species with the lowest number of GRs after cleaning and thinning was Guatteria anomala (with 85) and the one with the highest number was Bursera simaruba (with 2607). The number of GRs from each source (NFSI and GBIF) was variable for each species. Specific information is presented in Supplementary Material S2 (SM2). Information on bioclimatic variables and their selection, plus details concerning the six M areas used are included in Supplementary Material S3 (SM3). Warm-humid and warm-dry rainforest ecoregions were taken into account for species delimitation. Other ecoregions were considered when species (seven cases) were reported to be more widely distributed.

Niche Modeling and Species Distribution
For the species, the AUC ratio ranged from 0.730 to 0.952 and was considered acceptable (Supplementary Material SM5- Table S15). The best models of all species were significant (partial ROC = 0) and not very complex (delta AICc < 2). Only for T. amazonia and H. donnell-smithii were omission rates of 5.6 and 10.9%, respectively, obtained; the models of the other species had good predictive power (Error < 5%). In the final evaluation of the optimal models, there was excellent significance (partial ROC = 0) and variable omission rates (0 to 15.4%). The models obtained are better than a model given by chance (partial ROC = 0); their complexity and omission rates were variable. Applying GR thinning and discarding variables with the highest Pearson's correlation coefficient and less biological importance could be decisive in obtaining intermediate values of AUC (AUC = 0.5 is equivalent to a model given by chance), but a too high AUC (close to 1) is not an indication of a good model; in fact, it may be an indication of overfitting [25,52]. Applying the model used here resulted in models that perform better with new data if the model is good despite not having high AUCs [49,50]. The above is observed in the models obtained, because most did not present high AUCs but performed acceptably with new data (0 to 15.4% omission rates in the final evaluation phase (see Supplementary Material SM4)). The use of a single set of GRs may have been the reason for the highest omission value (15.4% for Vochysia guatemalensis). The training data set may be iterative (variant for one GR training group), for which replicates should be made with each GR group after which the average map obtained [63], which would have certainly improved the results. It has been found that the accuracy of the models is greater when considering species with small geographic ranges and limited environmental tolerance [62], that is, more accurate models are obtained from rare species than from widely distributed ones. Detailed information about the modeling is shown in Supplementary Material (SM4). For nine species, it was found that the variable with the greatest contribution was precipitation of wettest month (BIO13), for four species it was precipitation seasonality (BIO15) and for four other species it was temperature seasonality (BIO4) (Figure 1).

Mapping Potential Distribution and Changes
The current PD scenario of the species is shown in Figure 2. Dialium guianense has the smallest PD, while E. cyclocarpum and T. rosea have the broadest. The coincidence of the current PD area of our study with that obtained by Pennington and Sarukhán [1] was relatively low: the PD areas of M. mexicana, G. anomala and S. morototoni in our study differ with respect to the PD of the aforementioned authors by at least 50% of the area and have higher values, while the PD of G. glabra, P. campechiana and C. brasiliense of our study also show a percentage variation but have lower values. More details of the comparison can be found in Supplementary Material (SM5).
At least 14 species show the potential to occupy 97.4% of the territory currently oc cupied by the climax stage ( Figure 3). According to INEGI [4], the climax stage is distrib uted in 0.6% of the country (12,250 km 2 ). Other potential distribution areas are distributed in land use, mainly agricultural, pasture and different ecological succession of forests as shrub or herbaceous sprouts (Supplementary Material SM5- Table S16).

Mapping Potential Distribution and Changes
The current PD scenario of the species is shown in Figure 2. Dialium guianense has the smallest PD, while E. cyclocarpum and T. rosea have the broadest. The coincidence of the current PD area of our study with that obtained by Pennington and Sarukhán [1] was relatively low: the PD areas of M. mexicana, G. anomala and S. morototoni in our study differ with respect to the PD of the aforementioned authors by at least 50% of the area and have higher values, while the PD of G. glabra, P. campechiana and C. brasiliense of our study also show a percentage variation but have lower values. More details of the comparison can be found in Supplementary Material (SM5).
At least 14 species show the potential to occupy 97.4% of the territory currently occupied by the climax stage ( Figure 3). According to INEGI [4], the climax stage is distributed in 0.6% of the country (12,250 km 2 ). Other potential distribution areas are distributed in land use, mainly agricultural, pasture and different ecological succession of forests as shrub or herbaceous sprouts (Supplementary Material SM5- Table S16). Figure 4 shows the potential distribution of the 22 species, both for the assembly of all species (a) and for the three defined classes (b): Adequate-PD (14 to 22 species), Possible-PD (8 to 13 species) and Non-PD (equal or below 7 species). It can be seen that the Tuxtlas-Chimalapas region (southeast and center of the country, in the states of Veracruz, Oaxaca, Tabasco and Chiapas, number 3 and 6) is well represented by the assembly map. This is not the case for the Huastecas region (north of the state of Veracruz and south of the state of San Luis Potosí, number 1), nor for the Huasteca Potosi or Veracruzana region, as they are not represented in the expected class of the greatest representativeness.        Figure 4 shows the potential distribution of the 22 species, both for the assembly of all species (a) and for the three defined classes (b): Adequate-PD (14 to 22 species), Possible-PD (8 to 13 species) and Non-PD (equal or below 7 species). It can be seen that the Tuxtlas-Chimalapas region (southeast and center of the country, in the states of Veracruz, Oaxaca, Tabasco and Chiapas, number 3 and 6) is well represented by the assembly map. This is not the case for the Huastecas region (north of the state of Veracruz and south of the state of San Luis Potosí, number 1), nor for the Huasteca Potosi or Veracruzana region, as they are not represented in the expected class of the greatest representativeness. For the TR adequate-PD class, which represents a coincidence of at least 14 species, an occupation of 85,930 km 2 was found, equivalent to 4.45% of the country. According to INEGI [4], current land use is made up of 36% cultivated pasture, 29% at some successional TR stages, 15.5% annual rainfed agriculture, 1.9% water bodies, 1% urban and 17.6% other uses. The most considerable areas of intersection with other plant communities are with cloud forest in its shrub (2%) and old growth stages (1.8%) and with the pineoak forest arboreal stage (1.65%). For more details, see Supplementary Materials SM5.
Some differences in the PDs in our study with respect to those of Pennington and Sarukhán [1] may be due to the highly different performance by the Maxent and GARP algorithms [23,62]. Some of these differences are: the limit used in the definition of the binary layers (variable: determined by statistics; fixed: 0.3 presence probability), as well as the set of variables and GRs used in the creation of the models [40,42,62,64] and the periods in which it was carried out (2005 and 2019). However, there are coincidences towards the center of PD in most cases, possibly because the algorithms tend to be accentuated in the areas with a higher GR density, meaning that it is possible that the bias was not adequately reduced for all species, or that there are areas remaining to be sampled.
Climate change scenarios. The two climate change scenarios showed similarities in future PD areas ( Figure 5) and percent changes ( Figure 6). More information and details are presented in Supplementary Material SM5. For the TR adequate-PD class, which represents a coincidence of at least 14 species, an occupation of 85,930 km 2 was found, equivalent to 4.45% of the country. According to INEGI [4], current land use is made up of 36% cultivated pasture, 29% at some successional TR stages, 15.5% annual rainfed agriculture, 1.9% water bodies, 1% urban and 17.6% other uses. The most considerable areas of intersection with other plant communities are with cloud forest in its shrub (2%) and old growth stages (1.8%) and with the pine-oak forest arboreal stage (1.65%). For more details, see Supplementary Materials SM5.
Some differences in the PDs in our study with respect to those of Pennington and Sarukhán [1] may be due to the highly different performance by the Maxent and GARP algorithms [23,62]. Some of these differences are: the limit used in the definition of the binary layers (variable: determined by statistics; fixed: 0.3 presence probability), as well as the set of variables and GRs used in the creation of the models [40,42,62,64] and the periods in which it was carried out (2005 and 2019). However, there are coincidences towards the center of PD in most cases, possibly because the algorithms tend to be accentuated in the areas with a higher GR density, meaning that it is possible that the bias was not adequately reduced for all species, or that there are areas remaining to be sampled.
Climate change scenarios. The two climate change scenarios showed similarities in future PD areas ( Figure 5) and percent changes ( Figure 6). More information and details are presented in Supplementary Material SM5.    (loss refers to the area currently occupied by a given species that will no longer be suitable in the future; gain is for new areas outside the current distribution that will be suitable in the future). (loss refers to the area currently occupied by a given species that will no longer be suitable in the future; gain is for new areas outside the current distribution that will be suitable in the future). Figure 7 shows the projected changes for climate change assemblies for both scenarios. Additionally, the expected changes by number of matching species are shown (c and d). Considerable decreases in the TR Adequate-PD category are marked in the Tuxtlas-Chimalapas, Soconusco and Southern Tabasco regions, which are larger in these areas than the others (Figure 4). The increase in PD is mostly seen in the higher elevation regions, in the Huastecas and the Nautlas, while greater prevalence is indicated in the Lacandona rainforest region. The GFDL_CM3 model suggests higher humidity in the south of the country, while for the HADGEM2_ES model, they are more restrictive. However, those regions where more than 14 species coincide are considered optimal for conservation (dark green in c and d).
In the HADGEM2_ES scenario, decreases in the TR Adequate-PD category are projected to be 3.5% nationally. A similar situation is observed with the GFDL_CM3 scenario (see Supplementary Material SM5). According to the latter scenario, there are differences in areas of prevalence and loss of PDs, which are presented in Figure 8a. This means that two climate change scenarios find future coincidences in environmental variables that define potential distribution of 22 TR species. d). Considerable decreases in the TR Adequate-PD category are marked in the Tuxtlas-Chimalapas, Soconusco and Southern Tabasco regions, which are larger in these areas than the others (Figure 4). The increase in PD is mostly seen in the higher elevation regions, in the Huastecas and the Nautlas, while greater prevalence is indicated in the Lacandona rainforest region. The GFDL_CM3 model suggests higher humidity in the south of the country, while for the HADGEM2_ES model, they are more restrictive. However, those regions where more than 14 species coincide are considered optimal for conservation (dark green in c and d). In the HADGEM2_ES scenario, decreases in the TR Adequate-PD category are projected to be 3.5% nationally. A similar situation is observed with the GFDL_CM3 scenario (see Supplementary Material SM5). According to the latter scenario, there are differences in areas of prevalence and loss of PDs, which are presented in Figure 8a. This means that two climate change scenarios find future coincidences in environmental variables that define potential distribution of 22 TR species.

Proposed Conservation Areas
The conservation of 16,850 km 2 is proposed because the PD there is expected to be conserved for at least 14 species according to both climate change scenarios without considering the species emigration (Figure 8b). These are also the places where the best conserved TR relics are located in Mexico [3], that is, the Lacandona rainforest and the regions of greater elevation than the current PD area and where TR actually exists according to official cartography [4].

Niche Modeling and Species Distribution
This study considered the species with the highest structural index and those fre-

Proposed Conservation Areas
The conservation of 16,850 km 2 is proposed because the PD there is expected to be conserved for at least 14 species according to both climate change scenarios without considering the species emigration (Figure 8b). These are also the places where the best conserved TR relics are located in Mexico [3], that is, the Lacandona rainforest and the regions of greater elevation than the current PD area and where TR actually exists according to official cartography [4].

Niche Modeling and Species Distribution
This study considered the species with the highest structural index and those frequently mentioned in the literature. Such species resemble what in conservation biology are known as umbrella species (subset of the surrogate species), and their use in plant communities seems to be promising, since they are frequently used for modeling and protecting animal species. Despite being widely distributed [2,17,65], these species were considered in this study because they are a structural part of the plant community, and it would not be functional without some of these elements [66]. Furthermore, the transitions between successional stages occur gradually [57] and the presence of these species may be a reflection of this condition. Just a few species are used in planning the conservation of larger food webs to ensure that other coexisting species are also conserved. In fact, it is often used [67] in defining the minimum size of the area to be conserved, although this does not guarantee the conservation of rare or human interest species [68]. In this study, Ampelocera hottlei (Standl) Standl and Trophis racemosa (L.) were found to be important as they presented high rates of structural importance as in other TR studies [30].
A high percentage (85%) of the current area with TR fell into the Adequate-TR category. At least 14 species have the potential to occupy almost all of the current TR area (97.4%), which is enough to represent the old growth community and plan its conservation. It is possible to consider more species in further studies (e.g., in the Huasteca region), but the improvement in models and mapping may not be significant. Thus, we consider that our results contribute to TR conservation planning, particularly in old communities as well as in early ones.
The PD of the species is concentrated in the regions with the largest current surface area, that is, in the southern regions, called the Chimalapas and the Lacandona rainforest (see Figure 4b). However, there is a greater level of PD in the regions dedicated to agricultural and livestock use (52%), which coincides with the deforestation pattern experienced by TRs in Mexico during recent decades [2,69]. In regions where TR existed, several grass species were cultivated; however, there were also native species, which maintained the herbaceous successional state due to the continuous use of fire by its inhabitants to accelerate the re-growth of grasses [2]. This prevented the development of later successional stages [57]; this condition still persists [70]. If this alteration ceases, the community could repopulate the surrounding areas. The variants of this plant community are very diverse, as are those of the emerging secondary vegetation [2]; perhaps this is the reason why the current distribution area was not fully covered, which, although it was not this study's objective, was taken as a guide to measure the degree of certainty of the results.
The areas overlapping with other communities occur with the cloud forest shrub phase (2%) and with its old growth stage (1.8%), as well as with the pine-oak forest arboreal stage (1.6%). The overlaps can help in identifying ecotones, places where there is a combination of plant elements from different communities and which normally have greater floristic diversity [59,71]. The communities and their successional stages have no well-defined boundaries and they happen gradually; moreover, they are the result of their environment and the interactions between their biotic and abiotic elements.
While it is true that obtaining a more refined approximation of a species' PD requires considering interactions (including exclusionary interactions, not just symbiosis), our results approximate the original distribution because the methodology included: the interactions (biotic component, B); the environmental variables (abiotic component, A) and the M or accessible area (component M). Thus, the BAM diagram [47] was almost completely covered, indicating the conditions necessary for a species to live.

Projected Changes and Proposed Conservation Areas
The climate change scenarios consistently showed reductions in future PD. In small regions, there is an improvement in climate variables and, consequently, an increase in the conditions for TR growth. Moreover, the pressure and trend observed for species in humid tropical ecosystems is towards fragmentation and reduction. Hence, the scenario worsens for rare or small distribution species [7,72,73]. This pattern was reproduced both in individual species modeling and in the cartographic assembly of the species studied. Reduction in the largest current distribution areas and a possible shift to higher regions are expected. The PD areas in the TR Adequate-PD category (14 to 22 species) of the future scenarios were located in the areas that currently correspond to cloud forest (also in the Lacandona rainforest in the state of Chiapas). When compared to similar studies carried out for cloud forest [9], it is possible to foresee a more accentuated combination of floristic elements between both communities due to the change in their PD. However, in this study, we did not consider integrating the future adaptability of trees, appointed as the ability of forest species to grow successfully under climatic conditions different from those of their natural distributions [74].
In the base scenario and for most of the models, a significant contribution of the extreme values was found in the explanation of the model. Some cases referred to in the literature only refer to the TR from annual averages, that is, an average annual rainfall greater than 2000 mm with up to three or four months with 60 mm, an average annual temperature from 22 to 26 • C but not less than 20 • C and an average monthly temperature above 18 • C with fluctuations from 5 to 7 • C. They also refer to the "tropical rainy" climate that occurs in the TR [1,75]. However, our results indicate that intra-annual attributes are significant, such as the rainfall of the wettest month and its seasonality. We also found that the rainfall of the wettest month explains the PD of at least nine species studied. The above coincides with what Miranda and Hernandez [2] stated about TR distribution: not only is abundant precipitation necessary, but also its availability throughout the year. With increasing global warming and climate change, TR ecosystems may be at risk of abrupt and irreversible changes [11]. It is, therefore, important to have results that show the risk of exchange in conditions suitable for the TR. Many species will be unable to track suitable climates. As appointed by the Intergovernmental Panel on Climate Change -IPCC [11], the maximum speed at which a tree species can move across landscapes is below 20 km per decade.
Ecological niche studies are important to understanding the ecology of TR. Potential distribution shows potential to be used as complement research in conservation studies. The choice of proposed conservation areas took into account the long term since it is the most suitable for modeling changes in natural systems [12]. In addition, consensual PD areas were considered for two climate change scenarios. Maps of future PD show the likely distributions of TR and areas particularly vulnerable under climate change. Thus, areas where environmental conditions may prevail were located. Hence, if these areas are protected from now on, they can serve as resilience areas for biotic elements of the TR and, therefore, continue as the floristic refuge of countless species that have been able to survive until today. It should be considered that the processes of speciation, natural selection and the human race itself can have negative effects on some of its biotic elements, but here, TR conservation at the organizational level of plant community is sought. To achieve proper planning, one should not lose sight of and give priority to rare elements that are at some risk of disappearing from the natural environment. Studies must be accompanied by field verification to evaluate tree vigor and natural regeneration [76]. Managing the risks of climate change involves adaptation decisions with implications for future generations, as the promotion of new protected areas.

Conclusions
The potential distribution of species with ecological value in the tropical rainforest was modeled. In this process, significant correlations between geographic records and bioclimatic variables were found. Our results are consistent with the current distribution areas of the forest and suggest that with climate change scenarios it could be reduced. Therefore, it is possible to validate the proposed TR conversation areas in Mexico.
There are many challenges involved in modeling the distribution of a species, but there are even more challenges when modeling an ecosystem. The assembly we made allowed us to get closer to results for the tropical rainforest, but we are far from completing its modeling. As more accurate algorithms continue to be developed, the gap for modeling ecosystems and their related landscapes will narrow. With the currently available techniques and tools, we believe that this is not yet possible, since the rainforest is a diverse and complex community. Our work approaches using machine learning methods in an intensive way and rests on the pillar of having selected species with the highest structural index.
The challenges for forests are even greater since they are under high anthropic pressure. The fragmentation that this community is undergoing puts its future permanence at risk. In addition, more than half of the current potential area that we found is land being used for agriculture or livestock. For this reason, it is essential to conserve and manage it in a sustainable manner. Consideration should be given to where appropriate environmental conditions are expected to prevail when considering climate change scenarios.
We suggest that future research should pay close attention to the choice of species, explore the feasibility of modeling the community according to its successional stage and then assemble it, consider interactions between species and directly verify the results in the field.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to storage issues.