Predicting the Future Distribution of Ara rubrogenys , an Endemic Endangered Bird Species of the Andes, Taking Into Account Trophic Interactions

: Species distribution models (SDMs) are commonly used with climate only to predict animal distribution changes. This approach however neglects the evolution of other components of the niche, like food resource availability. SDMs are also commonly used with plants. This also suffers limitations, notably an inability to capture the fertilizing effect of the rising CO 2 concentration strengthening resilience to water stress. Alternatively, process-based dynamic vegetation models (DVMs) respond to CO 2 concentration. To test the impact of the plant modelling method to model plant resources of animals, we studied the distribution of a Bolivian macaw, assuming that, under future climate, DVMs produce more conservative results than SDMs. We modelled the bird with an SDM driven by climate. For the plant, we used SDMs or a DVM. Under future climates, the macaw SDM showed increased probabilities of presence over the area of distribution and connected range extensions. For plants, SDMs did not forecast overall response. By contrast, the DVM produced increases of productivity, occupancy and diversity, also towards higher altitudes. The results offered positive perspectives for the macaw, more optimistic with the DVM than with the SDMs, than initially assumed. Nevertheless, major common threats remain, challenging the short-term survival of the macaw.


Introduction
The problems of habitat destruction and climate change are the main threat to tropical mountain birds. Mountain bird species in the tropics are particularly at risk because they are isolated by hotter lowland zones which often makes them sedentary. In addition, when shifting their distribution up to higher altitudes, the new area of occupancy narrows [1]. The structure of the mountains itself also appears to be a constraining factor limiting the distribution shift with possible decline of habitat quality, for instance, the absence of suitable nesting sites or even vertical gaps between actual and potential future areas of distribution [2].
The slopes of the Andes are recognized as supporting the highest avian diversity in the world combined with high endemism rate but also more than 20% of threatened species [3]. In Bolivia, the red-fronted macaw (Ara rubrogenys Lafresnaye, 1847) is one of the 15 endemic species of this country [4]. Less than 30 years ago, A. rubrogenys was a little-known species [5]. It lives on the east Andean slopes of south-central Bolivia from 553 m up to 3094 m a.s.l. (Figure 1) and breeds between 1188 and 2696 m [6]. Its natural habitat is mainly semi-deciduous dry forest but this is most often severely degraded by pastoralism and by timber extraction into thorny scrubs with scattered trees [7]. While it was estimated that the threats were of limited extent in the early nineties [8], the status of the species worsens over the course of time with land conversion for agriculture, with poaching for illegal trade, with killing by farmers who consider them a pest and with poisoning with pesticides when they feed on the crops [6,[9][10][11]. The small breeding population (only 67 to 136 pairs) in eight close areas was also pointed out as major risk which increases their extinction risk due to correlated environmental fluctuations [6]. It was found that the birds use agriculture-scrub ecotones more than the forests for foraging, probably because the forests do not offer enough resources. A. rubrogenys is now ranked as "Critically Endangered" in Bolivia [11] and on the International Union for Conservation of Nature Red list [12]. In addition, it could be particularly threatened by climate change. In its area of distribution, A. rubrogenys uses only terrains along river valleys for roosting, feeding, resting and nesting [9]. Most of the nesting sites are located in steep river cliffs but such environments are not necessarily available at higher altitudes given the magnitude of warming predictions. Climate change is supposed to particularly affect the tropical Andes and notably Bolivia [13]; while warming already averaged 0.1 °C/decade between 1939 and 1998, it accelerated to 0.33 °C/decade between 1980 and 2005. Climate change scenarios suggest warming as high as 7.5 °C by 2080 and important modifications of the precipitation regime with respect to pre-industrial times. Species distribution models (SDMs) are based on the computation of an empirical relationship between the presence of a species (a sample of its distribution) and the actual values of the selected explaining factors. Range projections are obtained by computing a probability of presence over the study area using the relationship. The methodology has been applied to parrots several times with different objectives. For instance, comparing the projection of models fitted to historical data with actual ranges allowed to study the conservation status of Andean Pyrrhura in Colombia [14] or of Amazona in Venezuela [15]. Thanks to climate driven SDMs, the substitution of the actual climate with future climate allows projection of decreasing and shifting ranges in Amazona pretrei [16]. Climate driven SDMs also permitted computing of habitat suitability of 13 parrot species in invaded countries and to test the consequences of two successive trade bans in the US and the EU on the invasion success [17]. In some situations, however, climate factors are not sufficient explaining variables, they only improve the fit of models describing the niches in combination with other factors reflecting species requirements like habitat characteristics. Authors reached this conclusion for modelling Bonelli's eagle nesting sites with topography, disturbance, land-use or climate variable at several geographic scales [18]. Bad results were obtained with climate variables only, but climate significantly improves the quality of the prediction offered by the other sets of variables. In light of their results, the authors suggest that snow and low winter temperature may cause physiological stress hampering breeding success; however, they underlined that with more complex models, the interpretation of the effect of each explaining climate factor could become hard to find. Another interesting approach to model bird distribution consists of including biotic interactions as well as abiotic factors. For instance, in [19] SDMs driven by climate were used to simulate several shrub species making up the habitat of a bird species and the outputs of these models were set up as input variables to model animal presence with climate and topography. The authors found that this approach outperforms climate only models, which stressed the importance of taking into account, as far as possible, the different types of niche components to produce consistent simulations. This SDM approach was applied to refine the mapping of the suitability area of Amazona tucumana in Argentina and Bolivia. Here, the niche was defined with climate, land-use and the output of another SDM forecasting the distribution of a key plant resource for nestlings, Podocarpus parlatorei, providing niche cavities and food [20].
However, while SDMs driven by climate variables are now considered as a standard method to predict plant species distribution under future climate, this approach fails to consider the effect of the increasing CO2 concentration in air on plant physiology. Indeed, it is well established that increased air CO2 concentration improves the capacity of plants to resist water stress because the plants can minimize transpiration while still satisfying their CO2 requirements [21]. Contrary to SDMs, dynamic vegetation models (DVMs) are commonly able to reproduce this effect. Although, questions remain about the acclimation (organism trait responses occurring in days to weeks) and adaptation processes of plants (response occurring through evolutionary processes) to new climates which could lower future changes [22]. Forecasting under future climates with these types of models, forcing, or not, the increasing air CO2 concentration, gives contrasting results and the projected distributions of plant species or habitats under future climate with increased CO2 concentrations appear better conserved than with SDMs (e.g., [23,24]).
The objectives of this study are to evaluate the potential impact of climate change on the distribution of A. rubrogenys. We compare the results produced by SDMs driven by climate variables for A. rubrogenys and for 17 resource plant species over the area of A. rubrogenys and the results produced by a DVM for the same plant species, under present conditions and future climates (2070-2100) under RCP2.6 and RCP8.5 forcing. We predict that for the future, approaches with SDMs should conserve less of the original distribution area of A. rubrogenys than obtained with DVM.

Species Occurrences
Censuses of A. rubrogenys were conducted between 2008 and 2010 by two of us (E. Rocha Ledezma and L. Zúñiga Zeballos). The birds were observed between 5 and 11 a.m., in the course of linear surveys and near known nesting, feeding and roosting sites, while prospections were also realized in fields, forests and scrubs. Observations of A. rubrogenys in its area of distribution indicated that it mainly feeds on the fruits and seeds of the following wild species [25,26] (Table 1). A second check for duplicates was conducted after combining with climate factors for coordinates belonging to the same pixels. This had limited consequences on plant sample size but not for A. rubrogenys sample size which dropped to 63. A. rubrogenys' area of distribution was defined with the alpha-hull polygon method [27] using the ashape function of the R package "alphahull" [28].

Climate Data
For the current climate (monthly values for temperature, difference between maximum and minimum daily temperatures, precipitation, relative humidity, sunshine hours, wind speed), we used Worldclim version 2 for the time period between 1970 and 2000 at 2.5 arc-minutes [29]. For the future climate, we used the CMIP5 projections of the HadGEM2-AO global circulation model [30] under the representative concentration pathways of greenhouse gases corresponding to end of 21st century radiative forcings of 2.6 and 8.5 W/m 2 (RCP2.6 and RCP8.5, [31]). RCP2.6 would require a decline of greenhouse gas to reach no emission after 2072. With respect to the pre-industrial period, this scenario would keep global temperature rise below 2 °C in 2100. RCP8.5 is the worst hypothesis with emission of greenhouse gases continuing to increase throughout the 21st century producing global temperature increase between 2.6 and 4.8 °C for the period 2081-2100 [32].

Dynamic Vegetation Modelling
The DVM CARAIB (CARbon Assimilation In the Biosphere) was mainly described in [23,[33][34][35][36]. This model was initially conceived to simulate vegetation at global or continental scale and its response to climate change in the future or in the past [37][38][39][40][41][42][43]. The model was also applied to agricultural systems [44][45][46] or to tree species [24,[47][48][49]. It is a grid point model composed of 5 main interacting modules (hydrology budget, photosynthesis and stomatal regulation, carbon allocation and growth, heterotrophic respiration and carbon dynamics in the soil, competition between ecosystem strata and biogeography) and it is also possible to activate natural fire and migration modules. Input data are spatial monthly climates (minimal and maximal temperature, precipitation, relative humidity, sunshine hours and wind speed), CO2 air concentration, soil texture and color, elevation and a set of information describing the morpho-physiological characteristics of the plant species (traits), like the specific leaf area, leaf and sapwood C:N, plant height, deciduousness nature, et cetera, and climatic thresholds extracted from the distribution samples. Since trait information for the species was lacking, it was replaced by the values of the plant functional type to which the species belongs (Table 4 of Electronic Supplementary Material of [49]). Threshold values controlling germination and mortality under stress conditions are extracted from prescribed percentiles in their actual climate distribution extracted from the occurrence samples [36,40]. For computing the fitness statistics, we also ran the simulation for the present on the coordinates of the sets of pseudo-absences drawn for SDM modelling (see below). As output of the model, we used the net primary productivities of the species fractions (fNPP, gC/m 2 /y) and we selected thresholds of presence maximizing the true skill statistic (TSS, [50]). The threshold of presence allows us to compute the sensitivity (proportion of presences correctly predicted) and specificity (proportion of absences correctly predicted). We also computed the receiver operating characteristic curve (AUC). Positive TSS values and AUC larger than 0.7 indicate better agreement than random.

A. rubrogenys and Plant Species SDM
We used multiple logistic regressions, also called logit models. Studies have showed that it is difficult to rank the SDM methods according to their performances because they vary with dataset properties [51][52][53], but logistic regression has the best theoretical background and it is used here as a reference methodology for SDM. The models were computed in R with the glm function of the R package "stat" (The R Core Team) and we interpreted the output as a probability of presence (Px). Since absence data are also needed, we generated pseudo-absences for each species (A. rubrogenys and plant species). The pseudoabsences were sets of points randomly drawn around the occurrences containing twice as many points as occurrences. Our strategy is based on the fact that for rare events, rather than randomly sampling, it is more efficient to collect positive cases (presences), as much as possible, and to complete them with a limited number of random cases as absences [54]. There is no common rule to select the ratio of absences to presences and this varied considerably in ecological studies (see [55]). For each plant species, we built 10 datasets containing the occurrence coordinates and pseudo-absence coordinates drawn in radius buffers comprised of between 200 and 2000 km around presences. Furthermore, drawing pseudo-absences within the presences is justified by the fact that the species do not occupy their entire range. For each dataset, we identified the best combination of climate factor effects, in other words, a maximum of 6 (only 4 for N. herzogiana, due to smaller sample size) linear or quadratic effects, after exhaustive screening based on the Akaike information criterion (AIC), using the glmulti function of the R package "glmulti" [56]. This procedure allowed us to select the shortest distance for the pseudo-absences providing evident AUC increase, which may give more accuracy and meaningful fit of the models [57]. We finally used the datasets giving the highest AUC and the model selected in the above procedure with the lowest AIC. Possibly, we chose a model with slightly higher AIC so that most of the model coefficients may have Z-test p-values lower than 0.05. For validation we tested the selected models against their null models using the likelihood ratio test (with critical p-value = 0.05) and the AUC (auc function of "SDMTools" [58]). The cutting thresholds were also selected thanks to TSS. A linear effect alone reveals a strictly positive or negative ecological response to the considered factor. The combination of a linear effect with the quadratic effect of the same factor produces a bell-shaped ecological response curve. Owing to the logit link function, a quadratic effect alone acts as a threshold. As climate variables, we first selected annual mean temperature, temperature seasonality, maximum temperature of the warmest month, minimum temperature of the coldest month, annual precipitation, precipitation seasonality, precipitation of the wettest quarter and precipitation of the driest quarter. Then, to minimize collinearity, we computed the matrix of Pearson correlation coefficients in the species datasets and pointed the couple of variables with coefficients >0.7 [59]. Thus, we had to drop the maximum temperature of the warmest month, minimum temperature of the coldest month and annual precipitation for the plant species, and additionally, precipitation of the driest quarter for A. rubrogenys. For projections, we set missing value pixels with at least one of the explanatory factors outside the range encountered in the calibration datasets.
It should be noted that we tested the approach consisting of including outputs of the plant models to drive A. rubrogenys SDM using restricted estimation of the linear coefficients of the plant model output to impose positive coefficients. However, this approach failed mainly due to the fact that the productivity (DVM outputs) or the probability of presence (SDM outputs) were not uniformly higher over the A. rubrogenys area than outside. We understood this negative result as that animals are able to rely on resources made by plants growing in sub-optimal conditions, which seemed evident.

Results
Climate over the A. rubrogenys area varied a lot, exemplified by mean annual temperatures which fall between 12.31 and 23.00 °C and precipitation of the driest quarter which ranged between 6 and 62 mm (Table S1). Under RCP2.6 forcing, temperature and precipitation over the A. rubrogenys area increased while seasonality decreased with neverthe-less increases of the ranges between extreme higher and lower values. Under RCP8.5 forcing, temperature and precipitation followed the same direction but more intensively than under RCP2.6 forcing.

Dynamic Vegetation Modelling
Fitness statistics (Table 2) suggested limited to acceptable agreements for plant species DVM modelling. It is normal that the mechanistic model had significantly lower AUC than those obtained with the SDM. Indeed, first there is no training step in the DVM computation and, second, the AUC of the SDM has been maximized in the generation of pseudo-absences by varying the distance to the points of presence. Thus, comparing the performances of the two models in terms of AUC, or any other evaluator using pseudoabsences, cannot be made without bias. In this respect, it should be noted that the performance of the DVM in terms of sensitivity (Se) (which do not use pseudo-absences) is generally quite high, since they are generally larger than 0.8 or even 0.9. Otherwise, owing to excessive computing time, it was not possible to compute maps over the entire areas of distribution of the plant species. All the species but C. tubulosus were simulated as occurring over the A. rubrogenys presence area (Figures 2 and 3 and Figure S1). The range of fNPP also varied considerably between the species and the climate conditions. Table 2. Dynamic vegetation modelling (DVM) information: sample size for statistics (N), area under the receiver operating characteristic curve (AUC), true skill statistic (TSS), threshold for TSS (net primary productivity of the species fraction: fNPP, gC/m 2 /y), maximal net primary productivity of the species fraction over the area (Max fNPP, gC/m 2 /y), sensitivity (Se), specificity (Sp).   Figure S1).

SDM Modelling
It was possible to obtain models with acceptable to very good agreement for each of the plant species and for A. rubrogenys, with fitness indicators generally higher than those obtained for the DVM results (Table 3), except for the sensitivity (Se). For plants, projection maps revealed large areas of potential presence without any presence point which could be due to uneven sampling but also to natural barriers to migration, to competition or to speciation (Figures 4 and 5, Figures S2 and S3). Ten of the plant species were largely simulated as occurring over the A. rubrogenys' present area ( Figure 6). The simulated distributions over Bolivia sometimes showed large similarities with those produced by the DVM. For A. rubrogenys, the model defined a close region over the alpha hull defining the present area ( Figure 5). The original range was conserved under the future climates, the Px increased (Figure 6), while new suitable areas appeared connected with the former one ( Figure 5). Table 3. Species distribution models (SDM) information: sample size for model estimation and statistics (N), maximal distance from presence for pseudo-absences (D, km), area under the receiver operating characteristic curve (AUC), true skill statistic (TSS), selected threshold (Thresh.), sensitivity (Se), specificity (Sp). p-value of likelihood ratio test was lower than 10 −5 for each model; formula, i.e., polynomial parts of the logistic models are given in Table S2).

Comparisons of Plant Model Predictions for the Future
Under future conditions, changes predicted by the DVM over Bolivia were limited but plant presences tended rather to spread (Figures 2, 3 and S1). The trend was more pronounced over A. rubrogenys' area, particularly under RCP8.5 forcings while the fNPP of the pixel presence over the A. rubrogenys area behaved more or less similarly but with some exceptions (Cenchrus, P. visco, Figure 3). The mean number of plant species per pixel (plant diversity, Figure S4) computed for present climate was 10.3. It increased to 12.0 under RCP2.6 and 14.2 under RCP8.5 forcings. Furthermore, 80% (RCP2.6) to 96% (RCP8.5) of the pixels lost no species while a majority of them gained up to 11 new species (Figure 7). In addition, the maximal altitude increased for the majority of the species under the future climates (Figure 8).
With the SDMs, changes predicted over Bolivia were more significant than with the DVM (Figures 5, 6 and S3) but generally the presence also tended to spread. Over A. rubrogenys' area, no clear trend emerged for occupancy or Px. Diversity over pixels ( Figure  S4) for the present was higher than with the DVM (mean: 13.2) but decreased to 12.6 under RCP2.6 and to 13 under RCP8.5 forcings. Indeed, more pixels lost species while the gains were more modest (Figure 7). Furthermore, less species than with the SDM had maximal altitude increase under the future climates (Figure 8).

Discussion
The SDM allowed accurately modelling of A. rubrogenys presence with climate. Under future climate, occupancy, probability of presence and maximal altitude increased and new areas of occupancy appeared, connected with its present area. For the plants, DVM and SDM modelling allowed us to successfully compute the distribution of 17 plant species which are feeding resources of A. rubrogenys. Under present conditions, only one species, C. tubulosus, was not simulated as present over the area of A. rubrogenys by the DVM. SDMs gave better fits than the DVM, except for sensitivity (Se). The DVM mainly produced fNPP, and occupancy increases over the area of A. rubrogenys, also towards higher altitudes, under the RCP2.6 and RCP8.5 forcings. With DVM, we also obtained pixel diversity increases resulting from few species loss and more species gains. By contrast, the plant SDMs had no overall response under the future climate conditions for occupancy, probability of presence, altitude or diversity.
The lowest agreement of the DVM simulations compared to the plant SDMs were the result of lower specificity of the DVM. Here, we added an optimization step in which we selected a threshold for each species thanks to TSS. However, these thresholds were low compared to commonly fixed values. They probably reflected the fact that we used fNPP and thus the competition between species for water and light. Lowering the thresholds increases sensitivity at the expense of specificity. It was already noted in previous studies that the DVM may tend to simulate the fundamental niche of the species rather than the realized niche, in the absence of some critical biotic interactions. This could ultimately produce wider distributions compared to those produced by the SDMs as observed elsewhere [49,60]. Nevertheless, the results of both approaches, DVM and SDMs, were rather congruent for present conditions. For the future, the results were mostly in accordance with literature, in other words, the SDMs generally predict more widespread distribution shifts than DVM [24,49,[61][62][63] while this is not absolute. Therefore, it is not surprising that in limited parts of their present area, the Px of some species computed with SDM increased under the RCP forcings. In the Andes, SDMs predicted with dramatic reductions of the Andean vascular plant diversity for the future without excluding some small areas of better climatic stability where species numbers could increase [64]. Otherwise, the DVM results showing increases in diversity and the migrations towards higher altitudes seemed in accordance with the observations from the recent past. Permanent plots in the Andes showed thermophilization, in other words, the increasing number of tropical and subtropical tree species from lower altitudes showing shifts to higher altitudes [65]. This phenomenon, was also observed on European mountain summits. However, it should finally threaten the mountains species of the former communities because they are composed of slower growing species more adapted to harsher conditions compared to species from lower altitudes [66].
It is established that birds are directly sensitive to climate factors. While one of the main advantages of body temperature regulation in homeotherms is to allow an optimal functioning of the organisms under a wide temperature range, the limit of the mechanism is clearly the dissipation of heat in excess during exercise [67]. Thus, air temperature has to be considered as an effective component of the bird's niche. Homeotherms could also be directly sensitive to other climate factors than temperature, for instance, evaporation can limit the survival of birds in deserts during heat waves [68]. Furthermore, it was showed that birds have following their realized climate niche during the course of the last century [69]. For these reasons, the use of species distribution models (SDMs) with climate factors only to predict their future distribution has to be considered as a useful tool. The common response of birds in the Andes to climate change would be a decline [64,70]. Species with restricted distribution, like endemic species, are generally characterized by narrow climate niche and they seem to resist extinction by relatively high local abundance and good demographic resilience resulting from the accumulation of local adaptions [71]. The study of species, including parrots, naturalized outside their original area of distribution indicates however that the climate tolerance could be higher than estimated in the natural area, particularly for the species occupying narrow ranges of climate conditions or marginal climates in their native region [72][73][74]. With the climate factor, we obtained a very well-defined geographic range ( Figure 5). We supposed that the climate niche of A. rubrogenys would be broad compared to the other endemic species on the basis of the range of the climate factors over its present area provoked by the steepness of the climate gradient in the mountains (Table S1). Therefore, it would not be surprising that the SDM is driven by climate predicted occupancy increase of the present area under the RCP forcings and also on areas directly connected, at higher altitudes. However, this kind of situation would be rare. Another possibility would be that the area of A. rubrogenys is climatically stable enough to guarantee no species loss or even species gain even if those situations are rare [64]. If so, it would also explain the limited changes found for the plant species under the RCP forcings with both modelling approaches. Thus, as long as the nonclimate components of the niche of A. rubrogenys remain available, we could venture to forecast limited risk of extinction due to climate change provided that the diversity of the plant resources could increase and that species turnover could satisfy animal needs.
While the DVM approach is more promising than the SDM one for evaluating plant future, DVMs are nevertheless perfectible. DVMs need to elaborate an optimal strategy for parametrization [45,48,75]. It is indeed limited by the knowledge of the species-specific traits but also by the response of those traits to environmental factors and results could be improved by determining those responses and integrating them into the DVM [76][77][78]. The collection of plant traits with biogeographic information is a fast-growing field of knowledge [79]. Since the DVM simulates competition for light and water resources, another interesting point could be to integrate the species of the former communities with the potential newcomers to compute the emerging communities. The newcomers, in addition, could also constitute potential resources.

Conclusions
Under the milder RCP2.6 and the harsher RCP8.5 forcings, the SDM for A. rubrogenys and SDMs and the DVM for plant feeding resources of A. rubrogenys showed that the current area of suitability would be mostly preserved or even broadened. The actual feeding resources could be more precisely evaluated by taking into account the land-use and the productivity of the agricultural species which constitute important feeding resources for A. rubrogenys [6,25,26]. Nevertheless, predictions of the evolution of those resources in the course of time would require additional effort since uncertainties increase with the characteristics of the cultivated varieties and the decisions of the landowners or the farmers. The use of the DVM for predicting the future of plant species distribution and productivity is thus improvable. Our result boded well because A. rubrogenys would not be subjected to unprecedented climate conditions in the present area and the feeding resources so far would be preserved.
From topology, it could be possible to examine whether new suitable areas include cliffs appropriate for nesting. However, is it possible that the new suitable areas could be colonized even without appropriate cliffs thanks to behavioral flexibility? For instance, A. rubrogenys have been observed nesting in the palm Parajubaea torallyi [80] but, this tree species has an extremely small range distribution, and there are no other tree species that could offer large cavities for nesting within the range of the red-fronted macaw. Moreover, recent population genetic analyses show a low capacity for the species to colonize distant areas. Despite its restricted range, the population is structured in genetic clusters with low or null gene flow among them, despite that some colonies are separated by few tens of kilometers, that there are no ecological barriers and that macaws make larger daily and seasonal movements for feeding. Therefore, it is highly improbable that the species could colonize very distant areas in its future niche suitability [81].
Thus, while bird decline at global scale seems mainly provoked by increasing temperature [82], the main concerns about A. rubrogenys' future remained the direct human threats, in other words, habitat conversion, poaching and killing, which need to be urgently solved before the species rapidly goes extinct in nature [6,11,12].
Supplementary Materials: The following are available online at www.mdpi.com/1424-2818/13/2/94/s1, Table S1: Bioclimate variables over Ara rubrogenys area and changes under the RCP2.6 and RCP8.5 forcings, Table S2: Coefficients of the climate logistic models for Ara rubrogenys and the 17 plant species, Figure S1: Plant species distribution over Bolivia predicted by the DVM CARAIB for present and future under RCP2.6 andRCP8.5 forcings, Figure S2: Plant species continental distribution predicted by the SDMs for the present, Figure S3: Plant species distribution over predicted by the SDM for present and future under RCP2.6 andRCP8.5 forcings, Figure S4 Data Availability Statement: Information on data from publicly archived datasets used in this study is given in the Materials and Methods section.