Modelling Current and Future Potential Habitats for Plantations of Eucalyptus grandis Hill ex Maiden and E. dunnii Maiden in Uruguay

: Eucalyptus grandis and E. dunnii have high productive potential in the South of Brazil, Uruguay, and central Argentina. This is based on the similarity of the climate and soil of these areas, which form an eco-region called Campos. However, previous results show that these species have di ﬀ erences in their distribution caused by the prioritization of Uruguayan soils for forestry, explained by the particular conditions of each site. In this study, the site variables (climate, soil, and topography) that better explain the distribution of both species were identiﬁed, and prediction models of current and future distribution were adjusted for di ﬀ erent climate change scenarios (years 2050 and 2070). The distribution of E. grandis was associated with soil parameters, whereas for E. dunnii a greater e ﬀ ect of the climatic variables was observed. The ensemble biomod2 model was the most precise with regard to predicting the habitat for both species with respect to the simple models evaluated. For E. dunnii , the average values of the AUC, Kappa, and TSS index were 0.98, 0.88, and 0.77, respectively. For E. grandis , their values were 0.97, 0.86, and 0.80, respectively. In the projections of climatic change, the distribution of E. grandis occurrence remains practically unchanged, even in the scenarios of temperature increase. However, current distribution of E. dunnii shows high susceptibility in a scenario of increased temperature, to the point that most of the area currently planted may be at risk. Our results might be useful to political government and foresters for decision making in terms of future planted areas.


Introduction
In Uruguay, commercially, Eucalyptus trees have been planted in different areas of the country and today trees of this genus occupy 726,000 ha [1]. There is still a large area of soils prioritized for forestry by the government that could be planted, according to the current legislation (3.288 million ha, approximately). compared with simple models, and has been applied to project the distribution of species in different climate change scenarios [23]. Besides, some authors suggest the use of several types of algorithms (instead of just one) to predict the distribution of species since this allows to reduce uncertainty in marginal habitat situations [24]. Finally, SDM have been applied in this context of prediction of the potential distribution non-native species, particularly invasive non-native species [25]. Less frequent has been its use for introduced forest species with yield potential. However, those models may help us to understand the biological impacts and management practices for environmental suitability.
We hypothesized that the climate change that he expected climate change occurs in the region occupied by Uruguay has effect both on the area planted with the Eucalyptus species and on their productivity. The general objective of this work was to predict the current and future potential habitat of Eucalyptus grandis and E. dunnii in soils prioritized for forestry in Uruguay, under different future climate change scenarios. The specific objectives were (i) to identify the environmental variables that show the greatest association with the probability of occurrence of both species, (ii) to fit a prediction model of current and future habitat for the studied species, and (iii) to identify the regions with the greatest potential for the growth of these species.

Study Area
The study area corresponds to the plantations of Eucalyptus grandis and E. dunnii in the forest aptitude soils of the northern, western, and south-eastern regions of Uruguay ( Figure 1, Table S2 Supplementary Materials). Forest priority soils cover areas defined as forest according to current legislation, which include those marginal soils for agricultural uses [26]. use of several types of algorithms (instead of just one) to predict the distribution of species since this allows to reduce uncertainty in marginal habitat situations [24]. Finally, SDM have been applied in this context of prediction of the potential distribution non-native species, particularly invasive nonnative species [25]. Less frequent has been its use for introduced forest species with yield potential. However, those models may help us to understand the biological impacts and management practices for environmental suitability. We hypothesized that the climate change that he expected climate change occurs in the region occupied by Uruguay has effect both on the area planted with the Eucalyptus species and on their productivity. The general objective of this work was to predict the current and future potential habitat of Eucalyptus grandis and E. dunnii in soils prioritized for forestry in Uruguay, under different future climate change scenarios. The specific objectives were (i) to identify the environmental variables that show the greatest association with the probability of occurrence of both species, (ii) to fit a prediction model of current and future habitat for the studied species, and (iii) to identify the regions with the greatest potential for the growth of these species.

Study Area
The study area corresponds to the plantations of Eucalyptus grandis and E. dunnii in the forest aptitude soils of the northern, western, and south-eastern regions of Uruguay ( Figure 1, Table S2 Supplementary Materials). Forest priority soils cover areas defined as forest according to current legislation, which include those marginal soils for agricultural uses [26].  Its suitability is associated with the conditions which allow good forest growth: good conditions for root formation, adequate drainage, and low natural fertility. Currently, the total area of these soils is 4.42 million hectares and this land is characterized by having a low or average suitability for arable crops or livestock, according to the national productivity index.

Sources of Data and Environmental and Edaphic Variables
The data used in this analysis come from the National Forest Inventory developed by the General Directorate of Forestry for the years 2010, 2011, 2014, and 2016 [27], in which 128 parcels of Eucalyptus dunnii and 326 parcels of E. grandis were measured. The number of trees measured was: 5041, 1200, and 3561 in the first, second, and last inventory, respectively. The number of plots was determined to have a sampling error of less than 10%. The age range measured was from 5 to more than 20 years, the temporary plots were circular with an area of 113, 314, 616, and 1018 m 2 (6, 10, 14, and 18 m radius, respectively) and the number of trees in each plot ranged between 5 to 59 individuals. In each of the plots the diameter breast height as measured, the total and commercial height of all trees that exceed a height of 1.3 m and with diameters less than 10 cm and based on these data, survival was calculated. The diameter breast height was measured with a caliper with an accuracy of 1 mm, with two measurements perpendicular per tree and the height was measured with Vertex and Bitterlich's relascope [27]. These occurrence records were spatially filtered to avoid spatial autocorrelation [28] and spatial sampling bias [29,30]. There were generated ten sets of equal number, to occurrence for each species, of randomly distributed pseudo-absences [31] within the study area, with a minimum spatial distance of~1.5 km, larger than the hypotenuse of the pixel resolution, using the BIOMOD_FormatingData within the biomod2 R package [22]. Additionally, an equivalent random number of absence data were generated covering the area where the species are not present in Uruguay. Those areas were identified corresponding to plots of the National Forest Inventory without the presence of selected species.
A preliminary set of 19 bioclimatic and 48 monthly climatic variables were selected, either for the current situation or for future projections for the specific years 2050 and 2070. These were obtained from the Worldclim database [32] with a resolution of~1 km (Table 1). A model of global circulation of the atmosphere (GCM), which provides projections of the carbon dioxide concentration in the atmosphere [33] and considers a climate reconstruction (Community Climate System Model, CCSM4), and four representative pathway scenarios (RCP) of different greenhouse concentrations (2.6, 4.5, 6.0, 8.5) were used. We selected the CCSM4 climate change scenario because it is close to the ensemble average of whole GCMs both in terms of temperatures and rainfall and also because have been proof to present the highest accuracy to estimate regional temperature in north-eastern Argentina [34], close by our study area. The chosen values represent the increase in the heat absorbed by the Earth (for the year 2100) according to the concentration of greenhouse gases in each projection, measured in Watts per square meter. In addition, we used 20 edaphic (scale 1:40,000) and topographical variables [35], which were considered as constant for future projections.

Variable Selection
The environmental variables (n = 98) were reduced by the variable inflation factor (VIF < 10) [28] and the AUCRF R package [36,37]. We considered the most parsimonious model, i.e., the one with lower number of variables [28], built with the non-collinear variables that maximized the AUC value of the Random Forest (RF) model prediction [36,37]. This index was calculated as: where R 2 is the coefficient of determination. The selection of non-collinear variables was carried out using the stepwise procedure with the usdm package in the R program, depending on the importance of each variable [38] in the R program [39]. It was determined through simple linear correlations between the predictions of the model including all the variables (full model) and the prediction excluding the evaluated variable (reduced model) [28]. The AUCRF R package implements a stepwise variable selection procedure and returning the AUC value reached. The order in which the variables are added to the model is estimated according to the variable important measurement. The variable importance measurement offers two different importance measures, the mean decrease Gini (MDG) and mean decrease accuracy (MDA), respectively. Furthermore, the probability that the variable is elected to build the model. For more details see [36]. Hence, the model is firstly feed with the variable that returns the best AUC and new variables are added until the model is built with all available variables. The package returns the number and identity of the variables which give the best AUC value with their variable importance and probability of selection score [36].

Statistical Models
In this work, we used an ensemble approach (biomod2) [22] to describe the current and future habitat of Euacalyptus grandis and E. dunnii in Uruguay, as a basis to study potential changes in the optimal plantation areas for both species, under different climate scenarios. The biomod2 R package assembles the SDM which includes the Generalized Lineal Model (GLM), Generalized Additive Models (GAM), classification and regression trees (CART), Boosted Regression Tress (BRT), Random Forest (RF), Maximum entropy (MaxEnt) flexible determinant analysis (FDA), artificial neuronal networks (ANN), Multivariate Adaptive Regression Splines (MARS), and Surface Range Envelope (SRE). We used ensemble models calculated using the mean, median, coefficient of variation, confidence interval (inferior and superior), committee averaging (CA), and probability mean weight decay (WD) of the single model prediction. The CA was achieved by a binary (presence/absence) transformation using the threshold of the single model predictions. The threshold is the maximum score of the evaluation metric ("True Skill Statistics", TSS) for the evaluated dataset. Subsequently, the probability value of each pixel was calculated as the mean of the single pixel predictions. The WD ensemble modelling scaled the individual model predictions according to their accuracy statistic value (AUC) and the sum of all individual models [37,40]. We made ensemble predictions based on all single models' projections with an AUC > 0.90.

Selection and Validation of the Model
The evaluation and selection of the model correspond to the determination of the predictive capacity based on the quantification of the error and the misclassified data. These errors may be of commission, which is the classification as an absence of a data point that is present (false negative), and of omission, which is the opposite. We randomly split our dataset into two subsets, 70% of the data to train the model and 30% for model evaluation using cross-validation (100) which yielded 100 different fits per model. Models with higher mean values and smaller variations were considered as being the most accurate ones [28]. The criteria used to determine the predictive capacity of the model were the Kappa coefficient (K), the Area Under the ROC (Operational Characteristic of the Receptor) Curve (AUC) and the coefficient TSS. The first of these was used to estimate the veracity of the maps developed and is a qualitative measure of the concordance between the categoric predictors and describes the concordance rate between the observed and predicted values. The values of K vary from +1 to −1, where +1 represents a perfect classification and negative values represent a random fitment [41]. The AUC is the graphical representation of the commission errors on the horizontal axis (presences classified incorrectly, 1-specificity), and the omission errors (real presences which are omitted, sensitivity) are represented on the vertical axis for the whole value range [42]. The values obtained vary from 0.5 to 1, where 1 represents a perfect classification (presences and absences) and 0.5 represents a random prediction. According to Thuiller et al. [43], the fitted model may be classified as poor (AUC < 0.8), satisfactory (0.8 < AUC < 0.9), good (0.9 < AUC < 0.95), or very good (0.95 < AUC < 1). The coefficient TSS was applied to estimate the precision of the model with the presence/pseudo-absence data [44]. This compares the number of correct predictions except for those attributed to the random and, in the same way as the K index, considers the omission and commission errors. Its values oscillate between −1 and 1, where 1 indicates a perfect concordance and negative values indicate random behavior [42]. This index is defined according to the following expression: In this case, K index and TSS were calculated according to maximum TSS threshold [45].

Comparison of Current and Future Distributions
Current and future model predictions were obtained at a pixel resolution~1 km. Each pixel of the prediction contains the probability of occurrence of the species ranging from 0 (not probable) to 1 (highly probable). To better visualization of the probability of occurrence the probability values were classified in four categories (0-25%, low; 26-50%, moderately low; 51-75%, moderately high; 76-100%, high) for the present or future occurrence of both species in Uruguay. To calculate the variation in the area occupied in the future, compared to the area currently occupied by the two species, the probability maps were reclassified to 0 (absence) and 1 (present), using the same threshold applied to calculate TSS and Kappa. The total area of the potential distribution for each single projection was calculated by counting the number of pixels and multiplies it by the area of a pixel and the changed area was presented as percentage of total area loss (<100%) or gained (>100%) [30].

Selection of Variables
The results of the variable selection process revealed the importance of 4 and 5 variables out of the initial 98 variables for Eucalyptus dunni and for E. grandis model projection, respectively, as being non-collinear (n = 27; VIF < 10); and highlighted by the variable selection procedure run by the AUCRF R package. The selected variables were classified by importance (1 being the most important, and 5 the least important), considering the mean decrease Gini importance coefficient estimate with the AUCRF R package ( Table 2). The selected variables with greater predictive power for E. dunnii included: the depth of the A horizon, the highest and lowest temperatures of April and May, respectively, and the average temperature of the driest month. In the case of E. grandis, the most important variables were the percentage of clay, the depth of the A horizon, the isothermality, the percentage of silt, and the orientation.
The response curves analysis shows that the probability of occurrence for E. dunnii decreases as the temperature of the driest quarter increases and increases in soils where the A horizon is deeper (Figure 2A). The effect of the temperature of the driest trimester remains constant as it rises above 15 • C, while as the depth of the A horizon increases above 4 cm a constant increase in the probability of occurrence is observed. For E. grandis, the probability of occurrence is associated positively with the thickness of the A horizon and the aspect, but negatively with the percentage of clay in the A horizon ( Figure 2B). The probability of presence is highest when the proportion of clay is near 0% and remains constant when it exceeds 25%. Similarly, to E. dunnii, an A horizon thickness greater than 4 cm determines an important increase in the probability of occurrence. The orientation has a quadratic effect, with a higher probability of presence in the extremes of the range; the probability is highest for aspects over 300 • (north-west). On the other hand, the isothermality shows a polynomic relationship with the probability of occurrence, the probability being highest with levels close to 47.5%. The silt content in the A horizon has an effect very similar to that of the clay content. thickness of the A horizon and the aspect, but negatively with the percentage of clay in the A horizon ( Figure 2B). The probability of presence is highest when the proportion of clay is near 0% and remains constant when it exceeds 25%. Similarly, to E. dunnii, an A horizon thickness greater than 4 cm determines an important increase in the probability of occurrence. The orientation has a quadratic effect, with a higher probability of presence in the extremes of the range; the probability is highest for aspects over 300° (north-west). On the other hand, the isothermality shows a polynomic relationship with the probability of occurrence, the probability being highest with levels close to 47.5%. The silt content in the A horizon has an effect very similar to that of the clay content.

Model Selection and Validation
The predictive capacities of the models for the evaluated species, assessed through AUC, Kappa, and TSS, and their standard deviation are presented in Figure 3. Overall

Model Selection and Validation
The predictive capacities of the models for the evaluated species, assessed through AUC, Kappa, and TSS, and their standard deviation are presented in Figure 3.
are similar, though the accuracy of BRT for E. grandis were on the range of ANN and RF.
The predicted values obtained using ensemble models were higher than those given by individual models, for both species, except with the K index with the RF model (Table 3 and Figure  3). For E. dunnii, the average values of the AUC, Kappa, and TSS index were 0.98, 0.88, and 0.77, respectively. For E. grandis, their values were 0.97, 0.86, and 0.80, respectively. With both species there was a similar precision with the ensemble model, relative to individual models, whereas for E. grandis the precision was increased by using the ensemble model. Ensemble model overcame in accuracy single model predictions.  The predicted values obtained using ensemble models were higher than those given by individual models, for both species, except with the K index with the RF model (Table 3 and Figure 3). For E. dunnii, the average values of the AUC, Kappa, and TSS index were 0.98, 0.88, and 0.77, respectively. For E. grandis, their values were 0.97, 0.86, and 0.80, respectively. With both species there was a similar precision with the ensemble model, relative to individual models, whereas for E. grandis the precision was increased by using the ensemble model. Ensemble model overcame in accuracy single model predictions.

Current and Future Habitat Projection
The probability maps for the current habitat distribution show greater potential for both species in the northern and western areas ( Figure 4). Likewise, the south-west appears to be the area most suitable for Eucalyptus grandis, while for E. dunnii the probability of occurrence is greater in the west of the country. In terms of area, E. grandis is the species with the highest potential for occurrence values of these indices The prediction of the future occurrence of both species is shown in Figures 5 and 6. For E. dunnii, a drastic reduction in the species' habitat is predicted for the year 2050, with respect to the current situation, whereas for the year 2070 the reduction is less conspicuous. In both cases, the restriction of the occurrence of E. dunnii is greater in the scenarios in which greenhouse gases increase (RCP 8.5 vs. 2.6). This reduction would reach almost 100% with respect to the first of the mentioned scenarios (RCP 2.6), for all the studied regions. These tendencies are shown in Table 4 and Figure S1. A reduction of over 95% in the area of occurrence, even for the scenario of lowest temperature, was predicted by the models CA and WD for 2050. This reduction is almost total for the scenarios with higher greenhouse gas concentrations. The predictions for E. grandis show that the probability of occurrence will remain virtually constant for all the evaluated scenarios. On the other hand, the probability of occurrence for the two-time series does not show great changes for any of the four possible scenarios of temperature increment. These tendencies are shown in Table 4 and Figure S1. total for the scenarios with higher greenhouse gas concentrations. The predictions for E. grandis show that the probability of occurrence will remain virtually constant for all the evaluated scenarios. On the other hand, the probability of occurrence for the two-time series does not show great changes for any of the four possible scenarios of temperature increment. These tendencies are shown in Table 4 and Figure S1 (Supplementary Materials).

Variable Selection and Model Precision
The results of this study show that it is possible to establish a clear relationship between the probability of occurrence of both species and a reduced number of variables related to soil and climate parameters.
For E. grandis, the values of occurrence increase as a function of the percentage of clay content decrease and increase as the depth of the surface horizon increases. These soil texture and depth variables are closely related to eucalypts growth given their major influence on water availability [46,47], nutrients availability [48], and the volume of soil explored by the roots. The distribution of the eucalypt subgenus has been attributed to different responses to water availability. This, among other reasons, is due to the development capacity of a broad root system by modifying the ratio of root tissue to growth tissues and the depth of root penetration based on water availability [49]. For Uruguay, similar results have been reported [50], leading to the conclusion that the available water in the soil is one of the variables that better explain E. grandis growth. Likewise, [51] found a negative effect of the clay and silt presence on E. globulus growth. Isothermality, which indicates the temperature stability through the year [52], has also been reported as one of the variables explaining the growth of Eucalyptus cloeziana [53] and E. grandis [54]. The strong dependence of E. grandis on soil parameters was also verified by [55] evaluating plantations from 7 to 17 in the state of Paraná, southern Brazil.
The positive effect of north-west to north-east aspects on the probability of occurrence could be related to the greater exposure to solar radiation and the higher temperatures to which the planted trees in these orientations are exposed [56,57], although the positive effect of these variables arises from the interaction with variables such as soil moisture, temperature, or slope [58]. The relationship between orientation and radiation and the effects on the growth of eucalypts has been studied although according to Paton (1980) cited by [49] eucalypt in general are relatively less sensitive to the number of hours of light than to room temperature. However, [59] they argue that for the conditions of latitude such as those of Uruguay (30 • to 33 • S) in the north-northeast aspect there are greater hours of light which determines a greater period of photosynthesis as well as an occurrence of higher temperatures during the months of winter.
For E. dunnii the effects of thickness of the A horizon about the probability of presence can be explained by the same causes as those mentioned above. The positive effect of the temperature during April may be explained by the combined effect of a high relative water content in the soil (with greater precipitation than evapotranspiration) and intermediate temperatures that would promote tree growth. The increase in temperature during summer contributes to greater evapotranspiration, which may cause a deficit in the potentially available water in the soil during this period, although the precipitation has a relatively uniform distribution during the year [60]. This is consistent with the results obtained by [61] in the sense that the inclusion of climate variables with monthly records allows a greater predictive power of the species compared to the use of annual average values. According to these authors, this is due to the increasingly frequent occurrence of extreme weather events and therefore the importance of knowing the behavior of the species in such events.
In both species, higher values and less variation of the AUC were obtained with ensemble models compared to the use of individual models, which confirms the advantage of using the former to predict the occurrence of species [20,62,63]. Values of TSS above 0.85 (both with the assembled model and with the RF model) represent excellent predictive power [64]. These two models have shown very high capacities for predicting habitat in several species [37,65]. The work [66] obtained slightly lower values for this index (0.66 and 0.78) when analyzing the habitat prediction of E. sideroxylon and E. albens. According to the scale used by [67], several models in the group analyzed showed good to very good precision in habitat prediction, with the highest values corresponding to the RF model, for both species. This evaluation scale considers the following categories: <0.20 = poor, 0.21-0.40 = limited, 0.41-0.60 = moderate, 0.61-0.80 = good, and 0.81-1.00 = very good concordance. In general, greater values of this index were obtained for E. grandis, with all the evaluated models. ANN also returned acceptable accurate predictions with good-very good values for K, TSS and AUC, though it also presents a large variability on the predictions which might be linked to the hidden relationships built within the model [28] which suggest that the results must be interpreted with caution. Surprisingly, MaxEnt gave a moderate accurate model prediction when is one of the most use and accurate method to predict species distribution [68].

Current Potential Habitat and Future Projection
The results obtained confirm the hypothesis that the predicted climate change affects the probability of occurrence only for Eucalyptus dunnii since the occurrence of E. grandis remains unchanged in the long term. The areas with the greatest probability of occurrence of both species are associated with the deepest A horizon, with low amounts of clay and silt, corresponding to soils of the North of the country (Figure 4A,B). This kind of soil structure favors the availability of water (even when the rate of evapotranspiration is high, as in the summer), drainage, and a large exploration volume for the roots. The high probability of occurrence is also related to the high average temperatures of the area, which decrease from NW to SE by 4 to 5 • C. This favors growth during periods with greater potential water availability, such as fall [60]. Conditions that promote the growth of E. grandis combine deep, well drained, loamy to slightly sandy soil and an annual precipitation of 1000-1800 mm, in a temperature range from 8 to 36 • C, the optimum being 26 • C (Table S1 Supplementary Materials) [69,70]. This species has low tolerance of frosts and may grow at altitudes close to 900 m.a.s.l., although it has been cultivated successfully at higher levels. The natural distribution of the species (Australia) occupies a large area, with a latitude from 17 • to 36 • south, and covers different growing conditions. The areas with a greater probability of occurrence in Uruguay, according to our work, are similar to the ones described for E. dunnii ( Figure 4B). Soils with the required characteristics are also present, although over a smaller extension, in some central and western areas, and in some reduced areas in the north-west. The soils of south-western Uruguay show less potential because they are shallower-and therefore have lower water holding capacity and nutrients availability and have higher silt and clay contents. The relative heterogeneity of the probability of occurrence of this species that is depicted in the map (at a national level) is explained by the high heterogeneity of soils, resulting from the variety of geological parent materials. On the other hand, temperature increases are foreseen for this region, towards ranges more favorable for the growth of this species. The cardinal temperatures of E. grandis are reported to be 8, 25, and 36 • C (minimum, optimum, and maximum, respectively), while the average values in Uruguay are 12.9, 17.7, and 22.6 • C, respectively. Therefore, an increase in temperature could have a negative effect on growth by influencing the availability of water, but it is not expected that it will negatively influence the physiology of the species.
The stability of the area occupied by E. grandis can be explained by the fact that its presence is closely related to soil characteristics and the topography, and less so to climatic variables, such as the temperature ( Figure 6). For this reason, the average increase in temperatures projected for this region of 1 to 1.8 • C [9] would not imply changes in the species' area, although effects on its productivity may occur. The expected temperature increase (particularly in the north) would be accompanied by an increment in precipitation of 2.5% to 7%, which would have a positive effect on tree growth [17]. For these authors, the response of forest crops to pests and diseases would change in a climate change scenario [71], as well as for extreme climatic events [72] or an increase in the CO 2 level in the atmosphere [73].
The reduction in the area of occurrence of E. dunnii is predicted for the period 2000-2050 ( Figure 5); afterwards, the area would stay relatively unchanged for all the temperature increase scenarios (Table 4). Although one of the most important variables explaining the occurrence of this species is the thickness of the A horizon, the future presence seems to be heavily influenced by the temperature change. Despite the forecasted increases in temperature and average annual precipitation, increased variation between periods within years is also forecasted [9]. Such studies predict temperature increases during summer (1.2 to 1.8 • C) that will be greater in the period 2020-2050. This increment may explain the decrease in the presence of this species, given the negative effect of this variable on the probability of occurrence ( Figure 2). Temperature increases in summer could cause plant water deficiency, given the negative balance between evapotranspiration and precipitation. Despite these predictions, some authors have determined some capacity of eucalypt to adapt to situations of lower water availability through changes in the architecture of the root system [74]. These changes have been detected even in species considered as "less specialized" which could be a strategy of adaptation to future climate change. To this must be added the important genetic diversity of these species which determines a response potential that is observed in a wide geographical distribution [75].

Management Implications
SDM are important tools for forest management because these models can reliably predict current and future areas suitable for native and introduced species. This information is crucial in silviculture for the long-term timing of management options [76]. Furthermore, the species responses to environmental variables can provide information to help target for forest management programs at regional and national scales as has been proposed on this study.
The implications for future management are different for both species. With the species E. dunnii shows it is possible to visualize a greater potential in the northern and coastal zones in sites that simultaneously combine deep soils and comparatively higher temperatures in the country (Figure 2A). On the other hand, the results obtained with E. grandis shows a higher potential to show the advantages of plants to this species in relatively deep soils with little clay and silt content ( Figure 2B). This type of soil is concentrated in the coastal and north-north-east regions. Therefore, the effect of the orientation of the slope, although it shows the greatest potential in the north-north-east direction, lacks practical effects, since the plantation takes into account aspects of soil conservation such as contour lines. In turn, the probability values of occurrence of the species associated with this variable are of low magnitude to be taken into account.
Eucalypts plantations will continue to expand in Uruguay, growing in the coming years on soils for regions with forest potential (Figure 1), suitability soils that are still unexplored and this expansion is likely to be different for both species. For the projection of the future surface of E. dunnii, it must be taken into account that although the species comes from a relatively small area, it has managed to expand to a wide range of environmental conditions in substitution of E. grandis, [77]. However, the predictions of our model and other works carried out in Australia, China and South America show certain vulnerability of this species to changes in temperature, particularly increases in the average temperature above 19 • C (considered as the upper limit of climatic suitability) [77,78]. In this sense, future strategy should consider two aspects: (i) to have the widest possible genetic base (field trials and/or germplasm bank) so as to be able to select genotypes tolerant to the eventual increase in temperature, and (ii) have of interspecific hybrids with other species of eucalypts with greater ability to adapt to climate change. With respect to the projection of E. grandis, it is expected that it will continue covering an important area (higher than the current one) with the probable scenarios of temperature increase. The wide distribution of that this species shows worldwide, adding to the possibilities of genetic improvement in favor of more extreme weather conditions [8], determine in the long term an optimistic outlook for foresters. The model projections obtained in this study confirm the climatic adaptability of the species and therefore its relative stability against climate change.
Finally, the projections for both species should be taken with caution (particularly with E. dunnii) since the future increase of temperature explained by an increase in the concentration of atmospheric CO 2 could cause changes in the photosynthetic rate and in the use of water. Studies carried out on these topics with eucalypts species show contrasting results (see for example [8]).

Conclusions
The results obtained in this work indicate that the growth of E. grandis is associated basically with soil parameters, while that of E. dunnii shows a greater association with the temperatures of the fall and summer months. The direct relationship of the growth of both species with the depth of the surface horizon of the soil determines the importance of the choice of site with regard to obtaining high levels of growth. From this point of view, the potential of these species is greatest in the northern and coastal soils of the country, which, in general, have a comparatively greater volume of soil to be explored by the roots than the soils of the south-east of the country. With E. grandis, the positive effect of the north-north-west orientation on growth, but conditioned by soil conservation practices, must be considered. In this work, the ensemble models of habitat prediction emerge as useful tools to identify the most suitable areas for both species, based on the fitment values obtained. Eucalyptus grandis shows greater plasticity than E. dunnii regarding the different agroclimatic conditions of the country. Predictions of the future habitat of both species indicate that E. grandis is a species that could be used with certainty in the long term, in a wide variety of sites, whereas with E. dunnii, there may be areas of higher risk due to the probable climate change. Description of the potential current and future habitat allows selection of the areas with greater value for commercial plantations, which can assist afforestation plans, including the selection of appropriate genetic stock materials.
Despite the results obtained in this study, some limitations should be considered. Restrictions related to the anthropic distribution limits of the species are poorly known and could be insufficient to be effectively used in SDM. This is particularly important for the absence data, which may reduce the geographical extrapolations beyond the sampled area, leading to spurious results [29]. Thus, integrated and comparative survey strategies, including more detailed distribution information, should be considered in future studies. Additionally, an excessive simplification of the variables included in the model may mean that the selected variables limit their ecological interpretation (e.g., over-importance of the edaphic variables in our case). One possible way to avoid this problem is by making models that progressively include the set of environmental variables to interpret their importance step by step (see, for example [79]).