Species Ecological Envelopes under Climate Change Scenarios: A Case Study for the Main Two Wood-Production Forest Species in Portugal

: Species ecological envelope maps were obtained for the two main Portuguese wood-production species ( Eucalyptus globulus Labill. and Pinus pinaster Aiton) and projected future climate change scenarios. A machine learning approach was used to understand the most inﬂuential environmental variables that may explain current species distribution and productivity. Background and Objectives: The aims of the study were: (1) to map species potential suitability areas using ecological envelopes in the present and to project them in the future under climate change scenarios; (2) to map species current distributions; (3) to map species current productivity; and (4) to explore the most inﬂuential environmental variables on species current distribution and productivity. Materials and Methods: Climate, elevation data, and soil data sets were used to obtain present and future species ecological envelopes under two climate change scenarios. The o ﬃ cial land cover maps were used to map species distributions. Forest inventory data were used to map the species productivity by geostatistical techniques. A Bayesian machine learning approach, supported by species distributions and productivity data, was used to explore the most inﬂuential environmental variables on species distribution and productivity and to validate species ecological envelopes. Results: The species ecological envelope methodology was found to be robust. Species’ ecological envelopes showed a high potential for both species’ a ﬀ orestation. In the future, a decrease in the country’s area potentiality was forecasted for both species. The distribution of maritime pine was found to be mainly determined by precipitation-related variables, but the elevation and temperature-related variables were very important to di ﬀ erentiate species productivity. For eucalypts, species distribution was mainly explained by temperature-related variables, as well as the species productivity. Conclusions: These ﬁndings are key to support recommendations for future a ﬀ orestation and will bring value to policy-makers and environmental authorities in policy formulation under climate change scenarios.


WRBFU
Soil codes from the international soil classification system for naming soils and creating legends for soil maps The elevation data (Table 1) was derived from the SRTM (Shuttle Radar Topography Mission) with a spatial resolution of 30 m [31]. The soil data (Table 1) were derived from the ESDBv2 Raster Library in a grid of 1 km × 1 km [32][33][34].

Species Forests Cover (1995 and 2015)
To obtain maritime pine and eucalypts forests cover, the national reference thematic map for Land Cover and Land Use (LCLU) in Portugal named as COS (Carta de Ocupação do Solo) for 1995 (COS1995) and 2015 (COS2015) were used [35]. The COS maps have a scale of 1:25,000, one ha of the minimum cartographic unit, a five-level hierarchical classification with 89 (COS1995) and 48 (COS2015) classes in the most detailed level, available through Web Feature Services (WFS) [36], allowed a detailed analysis regarding these species' forests cover (Figure 1). To that end, the COS1995 legend was first harmonized concerning the COS2015 legend. In this analysis, only pure forests of the species were considered (i.e., forests where the dominant species represents more than 75% of the total number of trees and the ground cover is higher than 10%) as they represent more than 80% of these species' types of occupation [36,37].
S % Slope-The rate of change of elevation for each digital elevation model (DEM) cell (i.e., the first derivative of a DEM) A ° Aspect-The orientation of slope measured clockwise in degrees from 0 to 360, where 0 is north-facing, 90 is east-facing, 180 is south-facing, and 270 is west-facing.

WRBFU
Soil codes from the international soil classification system for naming soils and creating legends for soil maps The elevation data (Table 1) was derived from the SRTM (Shuttle Radar Topography Mission) with a spatial resolution of 30 m [31]. The soil data (Table 1) were derived from the ESDBv2 Raster Library in a grid of 1 km × 1 km [32][33][34].

Species Forests Cover (1995 and 2015)
To obtain maritime pine and eucalypts forests cover, the national reference thematic map for Land Cover and Land Use (LCLU) in Portugal named as COS (Carta de Ocupação do Solo) for 1995 (COS1995) and 2015 (COS2015) were used [35]. The COS maps have a scale of 1:25,000, one ha of the minimum cartographic unit, a five-level hierarchical classification with 89 (COS1995) and 48 (COS2015) classes in the most detailed level, available through Web Feature Services (WFS) [36], allowed a detailed analysis regarding these species' forests cover (Figure 1). To that end, the COS1995 legend was first harmonized concerning the COS2015 legend. In this analysis, only pure forests of the species were considered (i.e., forests where the dominant species represents more than 75% of the total number of trees and the ground cover is higher than 10%) as they represent more than 80% of these species' types of occupation [36,37].  Inventory (1990Inventory ( -1994 The National Forest Inventory data in 1990-1994 (NFI4) regarding the plots measured in pure stands (i.e., forests where the dominant species represents more than 75% of the total number of trees and the ground cover is higher than 10% [36,37]) of maritime pine and eucalypts were used to estimate Forests 2020, 11, 880 6 of 27 these species productivities ( Figure 2). The two-species stands characterization based on the NFI4 plots collected data [38] are synthesized in Table 2.  Inventory (1990Inventory ( -1994 The National Forest Inventory data in 1990-1994 (NFI4) regarding the plots measured in pure stands (i.e., forests where the dominant species represents more than 75% of the total number of trees and the ground cover is higher than 10% [36,37]) of maritime pine and eucalypts were used to estimate these species productivities ( Figure 2). The two-species stands characterization based on the NFI4 plots collected data [38] are synthesized in Table 2.  In Portugal, maritime pine stands are mainly established by natural regeneration and managed in medium to long revolutions (e.g., more than 35 years) [11,37,39]. Eucalypts stands are installed by plantation and commonly managed throughout three rotations of short harvesting cycles (e.g., 10-12-year coppices) [8,40].

Species Ecological Envelopes Maps
The distribution of maritime pine in Portugal (Figure 1a) corresponds to a coastal strip that extends from the Sado and the Tejo rivers basins to the Minho river basin, penetrating in the north and center inland up to elevations of 700-900 m, preferably in the southwest to north-facing slopes  In Portugal, maritime pine stands are mainly established by natural regeneration and managed in medium to long revolutions (e.g., more than 35 years) [11,37,39]. Eucalypts stands are installed by plantation and commonly managed throughout three rotations of short harvesting cycles (e.g., 10-12-year coppices) [8,40].

Species Ecological Envelopes Maps
The distribution of maritime pine in Portugal (Figure 1a) corresponds to a coastal strip that extends from the Sado and the Tejo rivers basins to the Minho river basin, penetrating in the north and center inland up to elevations of 700-900 m, preferably in the southwest to north-facing slopes where the Atlantic influence is still felt. Its climatic optimum corresponds broadly to an average temperature between 13 and 15 • C, an average temperature of the coldest month between 8 and 10 • C, and average annual precipitation between 500 mm and 1200-1400 mm. Regarding elevation, basal regions (0-400 m) are the most favorable. The species has severe limitations to growth above 800 m, due to wind and snow and has increased susceptibility to pests and diseases [40]. The species does not support well intense and prolonged colds and snow. As a pioneer species, it grows well in poor soil with light textures, preferably siliceous [39].
Eucalypts' present-day distribution in Portugal ( Figure 1b) includes areas of precipitation between the 700-2000 mm and mild winter due to its sensitivity to frost, in the north-coastal and central-coastal regions (Figure 1b) [38]. The species thrives in soil with a certain amount of humidity, well-drained, and without active limestone [7]. Data on eucalypts species naturalization in Portugal, showed that E. globulus is able to naturally regenerate from seeds both inside and outside plantations, corroborating that precipitation and distance to the sea (used as a surrogate of thermal amplitude) are key factors for the abundance of this species wildlings [9]. But topography, frost occurrence, and soil type also play a significant role [9]. Indeed, Catry et al. [9] observed that wildlings abundance, measured as the plant density through the number of plants per hectare, peaked at around 1500 mm of annual precipitation and it decreased with both low and high precipitation, reaching the lowest values below 800 mm and above 2400 mm. Higher plant density was also found in areas with milder temperatures, namely closer to the sea (with lower thermal amplitude) and with a lower number of days with frost. Plant density also seemed to be favoured in areas with intermediate elevation, higher slope and with certain soil types (namely, Cambisols and Podzols) [9]. These results are consistent to the previous knowledge for the species [38] and allow to understand what climate and site conditions this species natural regeneration prefers. Thus, helping to define the environmental variables range for the species site suitability. It should be noted that wildlings abundance is not a site productivity indicator and will be treated further in Section 2.2.3.
The species ecological envelopes for maritime pine and eucalypts [26], based on the species distribution information obtained in the literature [27], are defined by the following five variables: (i) temperature range (T max Aug-T min Jan); (ii) temperature limits (T max Aug and/or T min Jan); (iii) precipitation (P, i.e., BIO12); (iv) elevation (E); and (v) lithology. The variable thresholds (Table 3) ensure that more than 75% of the species forest inventory plots (IFN4), as a dominant species, are inside of the envelope. The maps of the input variables are reclassified by the Boolean method (0-Not suitable, 1-Suitable). The species ecological envelopes maps are obtained by map algebra, by summing the species variable maps produced earlier, resulting in five suitability classes classification as follows: (5) Excellent, (4) Good, (3) Regular, (2) Marginal, and (1) Unsuitable [25]. Legend: T max Aug-Maximum temperature in August ( • C); T min Jan-Minimum temperature in January ( • C); P-Annual precipitation (mm; BIO12); E-elevation (m).
In this study, the above-referenced methodology was used to produce the species ecological envelopes by using both the present climate data for 1960-1990 (Table 3 and Figure 3) and 1970-2000. The resulting envelopes were then applied using the predicted future climate data for 2050 and 2070 under the two RCPs scenarios (RCP 4.5-Meeting the Kyoto goal; and RCP 8.5-Not meeting the Kyoto goal). Afterwards, the ecological envelope classes area for both species in Portugal was extracted to analyze the country species potentiality in the present (1960-1990 and 1970-2000) and in future scenarios (2050 and 2070, both under RCP4.5 and RCP 8.5 scenarios).
Forests 2020, 11, x FOR PEER REVIEW 8 of 28 extracted to analyze the country species potentiality in the present (1960-1990 and 1970-2000) and in future scenarios (2050 and 2070, both under RCP4.5 and RCP 8.5 scenarios).

Species Forests Cover Distributions (1995 and 2015)
The main change trends regarding these two-species are well known. From 1995 to 2010, the maritime pine forest area decreased by 27% while eucalypts forest areas observed an increase of 13%, mainly by the reconversion of maritime pine areas after fire [11]. In this study, the land cover change map for the 20-year period (1995-2015) was performed, with a 1 km spatial resolution, by map algebra (subtraction) of the species forest cover in 1995 and 2015 (from COS1995 and COS2015).
The area of the ecological envelope classes over both species distributions in 1995 (COS1995) was extracted to analyze the potentiality of species forest areas (present and 2050 and 2070, both under RCP4.5 and RCP 8.5 scenarios). The ecological envelope classes change over both species distributions in 1995 (COS1995) was also performed by map algebra (subtraction) of the present and future scenarios (2050 and 2070, both under RCP4.5 and RCP 8.5 scenarios) to analyze the change trends in both species' potential areas.
The same grid (1 km × 1 km) was used to extract for each pixel the two-species percentage coverage in the COS1995 (Figure 1; Pb95 and Ec95; species presence occurs when pixel percentage coverage is higher than 0) and the correspondent bioclimatic , topographic and soil variables ( Table 1). The species presence was searched in every COS class where the species is mentioned. These data (Table A1 in Appendix A) were used in a machine learning approach to gauge variable influence in these two-species distributions ( Figure 4).

Species Forests Cover Distributions (1995 and 2015)
The main change trends regarding these two-species are well known. From 1995 to 2010, the maritime pine forest area decreased by 27% while eucalypts forest areas observed an increase of 13%, mainly by the reconversion of maritime pine areas after fire [11]. In this study, the land cover change map for the 20-year period (1995-2015) was performed, with a 1 km spatial resolution, by map algebra (subtraction) of the species forest cover in 1995 and 2015 (from COS1995 and COS2015).
The area of the ecological envelope classes over both species distributions in 1995 (COS1995) was extracted to analyze the potentiality of species forest areas (present and 2050 and 2070, both under RCP4.5 and RCP 8.5 scenarios). The ecological envelope classes change over both species distributions in 1995 (COS1995) was also performed by map algebra (subtraction) of the present and future scenarios (2050 and 2070, both under RCP4.5 and RCP 8.5 scenarios) to analyze the change trends in both species' potential areas.
The same grid (1 km × 1 km) was used to extract for each pixel the two-species percentage coverage in the COS1995 (Figure 1; Pb95 and Ec95; species presence occurs when pixel percentage coverage is higher than 0) and the correspondent bioclimatic , topographic and soil variables ( Table 1). The species presence was searched in every COS class where the species is mentioned. These data ( Table A1 in Appendix A) were used in a machine learning approach to gauge variable influence in these two-species distributions ( Figure 4).

Species Productivity Maps
In Portugal, the annual average maritime pine productivity (for fully mature trees of 35-45 years old) is about 5-10 m 3 ha −1 y −1 of timber in the north of the Tagus river, and lowers, in general, to 4 m 3 ha −1 y −1 in the south of the Tagus river. Much lower productivities of 0.7-3.4 m 3 ha −1 y −1 , are found in poor-quality areas, such as plots in coastal dunes [41]. In the coastal center regions maritime pine achieves the highest productivities (e.g., mean annual increment higher than 13.3 m 3 ha −1 y −1 , at the age of 30 years old), followed by the northern and central montane and submontane regions (e.g., around 10.3 m 3 ha −1 y −1 ) [42].
Eucalypts have higher productivities in the Atlantic seaboard (e.g., higher than 15 m 3 ha −1 y −1 ) followed by the Tagus valley region (e.g., between 10 and 15 m 3 ha −1 y −1 ; mean annual increment at the age of nine years old of eucalypts plantations without genetic improvement). In low productivity regions for this species, such as the inland of the country, productivities can be lower than 10 m 3 ha −1 y −1 [43].
To assess the two-species productivities, the most widely accepted site productivity indicator, the site index, as used [44]. As forest site productivity is very dependent on both environmental and biotic factors present at the site where the forest species is present [44], the most recent prediction

Species Productivity Maps
In Portugal, the annual average maritime pine productivity (for fully mature trees of 35-45 years old) is about 5-10 m 3 ha −1 y −1 of timber in the north of the Tagus river, and lowers, in general, to 4 m 3 ha −1 y −1 in the south of the Tagus river. Much lower productivities of 0.7-3.4 m 3 ha −1 y −1 , are found in poor-quality areas, such as plots in coastal dunes [41]. In the coastal center regions maritime pine achieves the highest productivities (e.g., mean annual increment higher than 13.3 m 3 ha −1 y −1 , at the age of 30 years old), followed by the northern and central montane and submontane regions (e.g., around 10.3 m 3 ha −1 y −1 ) [42].
Eucalypts have higher productivities in the Atlantic seaboard (e.g., higher than 15 m 3 ha −1 y −1 ) followed by the Tagus valley region (e.g., between 10 and 15 m 3 ha −1 y −1 ; mean annual increment at the age of nine years old of eucalypts plantations without genetic improvement). In low productivity regions for this species, such as the inland of the country, productivities can be lower than 10 m 3 ha −1 y −1 [43].
To assess the two-species productivities, the most widely accepted site productivity indicator, the site index, as used [44]. As forest site productivity is very dependent on both environmental and biotic factors present at the site where the forest species is present [44], the most recent prediction models including environmental parameters, for each species site index evaluation, were selected [45][46][47]. The site index assessment for the NFI4 maritime pine plots was performed using the empirical model developed by Nunes et al. [45], a multiple asymptotes polymorphic equation from Cieszewski [48], modified to include parameters related to climate and soil characteristics (P or BIO12, T or BIO1, WINTER or BIO6, and ST 1 , ST 2 , ST 3 ; Tables 4 and 5). The parameters were obtained from the climate data set from the WorldClim , the elevation data from STRM (Shuttle Radar Topography Mission), and the soil data from the ESDBv2 Raster Library (Table 1). Table 4. Site index models for the species maritime pine [45] and eucalypts [46,47].

Species Site Index Model
Maritime pine SI 50 -Dominant height at the reference age of 50 years  Legend: SI 50 -hdom (m) at the age of 50 years; SI 10 -hdom (m) at the age of 10 years.
The site index assessment for the NFI4 eucalypts plots was performed using the empirical model developed by Tomé et al. [46,47] that includes one parameter related to climate (pd-Number of precipitation days per year; Tables 4 and 5). Regarding the climate parameter pd (Table 4) the observations from Portuguese meteorological stations in three temporal 30-year series (1931-1960, 1951-1980 and 1971-2000) were used [49][50][51][52][53][54]. After, the pd climatic parameter was interpolated using the Inverse Distance Weighted (IDW) algorithm. The allocated weight was the local correspondent elevation (E) of the species NFI4 plots.
Afterwards, to obtain the species productivities maps (Maritime pine site index-SI 50 ; and Eucalypts site index-SI 10 ), the data spatial characterization was performed through a 2-step approach.
The variogram is a vector function used to compute the spatial variability of regionalized variables defined by the following equation: 2 (1) are the numerical values of the observed variable at points x i , and x i + h. The number of forming pairs for h distance is N (h). Thus, it is the average value of the square of the differences between all pairs of points in the geometric field spaced at h distance. The graphic study of the variograms obtained provides an overview of the spatial structure of the variable. The nugget effect (C 0 ) shows the behavior at the origin, whereas the sill (C 1 ) and the amplitude (a) define the variance used in the subsequent interpolation steps and the spatial influence of the variable, respectively.
Secondly, a Sequential Gaussian Simulation (SGS) was used as a stochastic approach used to map the two-species productivities' spatial distribution, through the construction of 100 scenarios. SGS starts by defining the univariate distribution of values followed by a normal score transform of the original values to a standard Gaussian distribution. Simulation of the obtained normal scores at grid node locations is done sequentially with simple kriging (SK) using the normal score data and a zero mean [57]. Once all normal scores have been simulated, they are back-transformed to original grade values. The result of a simulation process is a twisted version of an estimation process, which reproduces the statistics of the known data, making a realistic view, but supplying a low prediction behavior. Nevertheless, if multiple sequences of simulations are designed, it is possible to obtain more reliable probabilistic maps. The average image was computed, together with the associated standard deviation for the spatial uncertainty visualization [58].
Finally, the resulting productivity mean images (MI) and the spatial uncertainty representation, through standard deviation (SU), with a 2.45 km × 2.45 km resolution, were reclassified in five productivity classes using the statistical criteria of Jenks natural breaks [59].
Afterwards, for each maritime pine and eucalypts NFI4 plots (1990)(1991)(1992)(1993)(1994), the simulated site indices (SI 50 and SI 10 , respectively) and the correspondent climatic , topographic and soil variables were extracted (Table 1). These data (Table A1) were used in a machine learning approach to gauge variable influence in these two-species productivities ( Figure 4).

Variables Influence Analysis in Species Distribution and Productivity
Bayesian networks have become extremely popular for gaining deep insight into manifold problem domains [60]. Their inherently graphical structure built upon machine learning algorithms based on the Bayesian Network formalism leads to robust models optimal for exploring and explaining complex problems. In this study, four supervised Bayesian networks machine-learned from the data (i.e., presence-Pb95 and Ec95 and productivity-SI 50 and SI 10 ) were created to assess the most influential variables conditioning maritime pine and eucalypts distribution and productivity. To measure the variable influence regarding the target node (X) of each network, first, the relative mutual information (RMI) was calculated to quantify the percentage of information gained on X by observing Y where MI(X, Y) is the mutual information [61] in absolute terms and H(X) the entropy associated with the probability distribution of X. The entropy of a variable can also be understood as the sum of the expected log-losses of its states. Both metrics are expressed in bits and can be calculated as follows where H (X|Y) is the conditioned entropy associated with the probability of X given Y (how likely is X if you learn Y) and p(x) is the probability of x.
From the RMI analysis, a ranking was obtained regarding the significance of each variable to the knowledge of the target node. However, due to the high dispersion, the large number of variables might be introduced in the model, a second measurement was conducted with the most influential variables identified by the RMI analysis. The Bayesian models were machine-learned again in a simplified scheme and a contribution analysis was computed, showing the contribution of each variable towards the predicted mean of the target variable (i.e., presence-Pb95 and Ec95 and productivity-SI 50 and SI 10 ). The contribution of each variable is estimated using counterfactuals based on a specific neutral state that represents the neutral condition of the system. For modeling and graphical representation, the Artificial Intelligence software BayesiaLab 9 [62] was used.
The Bayesian machine learning analysis for the data (i.e., presence-Pb95 and Ec95 and productivity-SI 50 and SI 10 ) was carried out with the Augmented Naïve Bayes algorithm. The architecture of this algorithm is rooted in a naive architecture (the target node is the parent of all the other nodes), enhanced by the relations of the child nodes [63]. Cross-validation was used to estimate the prediction accuracy of the supervised machine learning models created. The relative significance (RS), as the ratio between the mutual information brought by each variable and the maximum mutual information, which corresponds with the most informative variable, was used to interpret these results. Due to the informative dispersion by the large number of variables introduced in the model, a closer look into the variables with the highest RS was done, which were taken to a second step to build a simplified Bayesian network (RS > 70% and RMI sum > 80%) In the end, the open-source software WEKA [64] was used to fit Machine Learning (ML) regression tree models by the algorithm Random Forest [65,66] to compare the quality of ecological envelopes variables as regressors to the ones selected by the Bayesian network analysis.

Species Ecological Envelopes under the Climate Scenarios
The ecological envelope maps analysis for both species in Portugal (Figure 5a,d) indicates that, in the present , the country's area potentiality for maritime pine is distributed as follows: 34% excellent, 42% good, and 22% regular. For eucalypts, the country's area potentiality is distributed where 31% is excellent, 36% is good, and 28% is regular. Therefore, in Portugal, there is a high potential for both species' afforestation (an area over 95%). where H (X|Y) is the conditioned entropy associated with the probability of X given Y (how likely is X if you learn Y) and p(x) is the probability of x. From the RMI analysis, a ranking was obtained regarding the significance of each variable to the knowledge of the target node. However, due to the high dispersion, the large number of variables might be introduced in the model, a second measurement was conducted with the most influential variables identified by the RMI analysis. The Bayesian models were machine-learned again in a simplified scheme and a contribution analysis was computed, showing the contribution of each variable towards the predicted mean of the target variable (i.e., presence-Pb95 and Ec95 and productivity-SI50 and SI10). The contribution of each variable is estimated using counterfactuals based on a specific neutral state that represents the neutral condition of the system. For modeling and graphical representation, the Artificial Intelligence software BayesiaLab 9 [62] was used.
The Bayesian machine learning analysis for the data (i.e., presence-Pb95 and Ec95 and productivity-SI50 and SI10) was carried out with the Augmented Naïve Bayes algorithm. The architecture of this algorithm is rooted in a naive architecture (the target node is the parent of all the other nodes), enhanced by the relations of the child nodes [63]. Cross-validation was used to estimate the prediction accuracy of the supervised machine learning models created. The relative significance (RS), as the ratio between the mutual information brought by each variable and the maximum mutual information, which corresponds with the most informative variable, was used to interpret these results. Due to the informative dispersion by the large number of variables introduced in the model, a closer look into the variables with the highest RS was done, which were taken to a second step to build a simplified Bayesian network (RS > 70% and RMI sum > 80%) In the end, the open-source software WEKA [64] was used to fit Machine Learning (ML) regression tree models by the algorithm Random Forest [65,66] to compare the quality of ecological envelopes variables as regressors to the ones selected by the Bayesian network analysis.

Species Ecological Envelopes under the Climate Scenarios
The ecological envelope maps analysis for both species in Portugal (Figure 5a,d) indicates that, in the present , the country's area potentiality for maritime pine is distributed as follows: 34% excellent, 42% good, and 22% regular. For eucalypts, the country's area potentiality is distributed where 31% is excellent, 36% is good, and 28% is regular. Therefore, in Portugal, there is a high potential for both species' afforestation (an area over 95%).
Likewise, when analyzing the ecological envelope maps in the present, but now with the climate data sets from 1970 to 2000, little change was observed and the high potentiality for both species' afforestation maintains.  Figure 5d), showed a decrease in the country's area potentiality (excellent, good, and regular areas) for both species. Maritime pine potential area is expected to decrease to 63-54%, respectively, in the RCP4.5 scenario and the RCP8.5 scenario. However, for eucalypts the decrease will not be so severe, ranging to 73-62%, respectively, in each scenario.
The ecological envelope maps for the projected future scenarios in 2070 are very similar to the ones in 2050. However, the maritime pine potential area is expected to decrease further to 59-51%, respectively, in the RCP4.5 scenario and RCP8.5. Similarly, for eucalypts, the decrease is expected to be 68-59%, respectively, in each scenario ( Figure A1).

Species Forests Cover Distributions Maps
Comparing the ecological envelope classes (Figure 5a) over the present species distributions in 1995 (COS1995; Figure 1) it was verified that 99% of maritime pine forest and eucalypts forest were found inside excellent, good and regular ecological envelope classes (Figure 6a). In fact, 62% of the  Figure 5d), showed a decrease in the country's area potentiality (excellent, good, and regular areas) for both species. Maritime pine potential area is expected to decrease to 63-54%, respectively, in the RCP4.5 scenario and the RCP8.5 scenario. However, for eucalypts the decrease will not be so severe, ranging to 73-62%, respectively, in each scenario.
The ecological envelope maps for the projected future scenarios in 2070 are very similar to the ones in 2050. However, the maritime pine potential area is expected to decrease further to 59-51%, respectively, in the RCP4.5 scenario and RCP8.5. Similarly, for eucalypts, the decrease is expected to be 68-59%, respectively, in each scenario ( Figure A1).

Species Forests Cover Distributions Maps
Comparing the ecological envelope classes (Figure 5a) over the present species distributions in 1995 (COS1995; Figure 1) it was verified that 99% of maritime pine forest and eucalypts forest were found inside excellent, good and regular ecological envelope classes (Figure 6a). In fact, 62% of the maritime pine forest is inside the excellent class (Figure 6a) compared to 52% of the eucalypts forest (Figure 6a). Moreover, according to the known land cover trend between these two species, the maritime pine forest has experienced a reduction between 1995 and 2015 which was mostly replaced by eucalypts new afforestation (Figure 6b).
Little change was observed for the ecological envelope maps in the present, obtained using the climate data sets from 1970 to 2000 instead (Figure 7a).
In 2050 future scenarios (Figure 5b,c), the maritime pine forest area (Figure 1a) will lose potentiality by decreasing from 99% to 84-72%, respectively, in each scenario (Figure 7b-d). Nevertheless, for the eucalypts forest area (Figure 1b) the potentiality loss will be lesser than for maritime pine, decreasing only to 91-85%, respectively, in each scenario (Figure 7b-d). Indeed, maritime pine forest in excellent areas class will experience a great loss in potentiality, from 62% to 19-11%, respectively, in each scenario (Figure 5b-d). Conversely, the eucalypts forest showed a much lesser decrease, from 52% to 32-22%, respectively, in each scenario (Figure 5b-d). The loss of maritime pine forest area potential is more severe in the north inland of the country (three-class dropdown) followed by the central north (two-class drop drown; Figure 7b-d). However, the loss of the eucalypts' forest area potential will occur in the inland (three-class drop drown) followed by the center (two-class dropdown; Figure 7b-d).
Once more, the projected future scenarios in 2070 are very similar to the ones in 2050. However, the maritime pine forest area (Figure 1a) will continue to lose potentiality by decreasing to 78-65%, respectively, in each scenario ( Figure A2). Again, the eucalypts forest area (Figure 1b) potentiality loss will be lesser than for maritime pine, decreasing only to 89-82%, respectively, in each scenario ( Figure A2). Maritime pine forest in excellent areas class is expected to lose potentiality to 15-6%, Moreover, according to the known land cover trend between these two species, the maritime pine forest has experienced a reduction between 1995 and 2015 which was mostly replaced by eucalypts new afforestation (Figure 6b).
Little change was observed for the ecological envelope maps in the present, obtained using the climate data sets from 1970 to 2000 instead (Figure 7a).
In 2050 future scenarios (Figure 5b,c), the maritime pine forest area (Figure 1a) will lose potentiality by decreasing from 99% to 84-72%, respectively, in each scenario (Figure 7b-d). Nevertheless, for the eucalypts forest area (Figure 1b) the potentiality loss will be lesser than for maritime pine, decreasing only to 91-85%, respectively, in each scenario (Figure 7b-d). Indeed, maritime pine forest in excellent areas class will experience a great loss in potentiality, from 62% to 19-11%, respectively, in each scenario (Figure 5b-d). Conversely, the eucalypts forest showed a much lesser decrease, from 52% to 32-22%, respectively, in each scenario (Figure 5b-d). The loss of maritime pine forest area potential is more severe in the north inland of the country (three-class dropdown) followed by the central north (two-class drop drown; Figure 7b-d). However, the loss of the eucalypts' forest area potential will occur in the inland (three-class drop drown) followed by the center (two-class dropdown; Figure 7b-d).
Forests 2020, 11, x FOR PEER REVIEW 15 of 28 respectively, in each scenario ( Figure A1). Eucalypts forest will continue showing a much lesser decrease to 28-13%, respectively, in each scenario ( Figure A1). The same analysis was performed regarding the species distributions in 2015 (COS 2015; Figure  1) but only a small percentage variation (ranging ±3%) was observed.

Eucalypts-Present change (1995-2015)
Future change (1995-2050)  Once more, the projected future scenarios in 2070 are very similar to the ones in 2050. However, the maritime pine forest area (Figure 1a) will continue to lose potentiality by decreasing to 78-65%, respectively, in each scenario ( Figure A2). Again, the eucalypts forest area (Figure 1b) potentiality loss will be lesser than for maritime pine, decreasing only to 89-82%, respectively, in each scenario ( Figure A2). Maritime pine forest in excellent areas class is expected to lose potentiality to 15-6%, respectively, in each scenario ( Figure A1). Eucalypts forest will continue showing a much lesser decrease to 28-13%, respectively, in each scenario ( Figure A1).
The same analysis was performed regarding the species distributions in 2015 (COS 2015; Figure 1) but only a small percentage variation (ranging ±3%) was observed.

Current Climate Species Productivities Maps
The two-species current productivities maps (Figure 8), obtained from the National Forest Inventory data in 1990-1994 (NFI4), show high productivity areas for maritime pine (SI 50 > 21.15 m) mainly in the west-northern and the center of the country (Figure 8a), and for eucalypts (SI 10 > 21.27 m) mainly in a western strip in the center and north of the country (Figure 8b). Inland regions are associated with high spatial uncertainty, while coast regions are associated with low spatial uncertainty, as expected due to those species' present distributions (Figure 1). The productivity classes, obtained by using the statistical criteria of Jenks natural breaks, are less restrictive than the ones found in the literature ( Table 5). The species highest productivities were found overlapping areas from regular to excellent ecological envelope classes, namely, 68% for maritime pine and 56% for eucalypts (Figures 6 and 8).

Current Species Productivities Maps
The two-species current productivities maps (Figure 8), obtained from the National Forest Inventory data in 1990-1994 (NFI4), show high productivity areas for maritime pine (SI50 > 21.15 m) mainly in the west-northern and the center of the country (Figure 8a), and for eucalypts (SI10 > 21.27 m) mainly in a western strip in the center and north of the country (Figure 8b). Inland regions are associated with high spatial uncertainty, while coast regions are associated with low spatial uncertainty, as expected due to those species' present distributions (Figure 1). The productivity classes, obtained by using the statistical criteria of Jenks natural breaks, are less restrictive than the ones found in the literature ( Table 5). The species highest productivities were found overlapping areas from regular to excellent ecological envelope classes, namely, 68% for maritime pine and 56% for eucalypts (Figures 6 and 8).

Variables Influence Analysis in Current Species Distribution and Productivity
The Bayesian machine learning analysis was performed for the two-species current distributions (Pb95 and Ec95) and productivities (SI50 and SI10). The cross-validation indicated that Bayesian models for species presence showed a higher quality of the problem representation than those for

Variables Influence Analysis in Current Species Distribution and Productivity
The Bayesian machine learning analysis was performed for the two-species current distributions (Pb95 and Ec95) and productivities (SI 50 and SI 10 ). The cross-validation indicated that Bayesian models for species presence showed a higher quality of the problem representation than those for productivity. The contingency table fit (CTF) values for Pb95 and Ec95 (above 85%) indicate that the joint probability distribution by the created networks is good. In turn, the CTF values for SI 50 and SI 10 are slightly lower, ranging from 70% to 80%. This circumstance can be appreciated when comparing the RMI for the species presence and productivity (Table A2). The most contributive variables (RS > 70% and RMI sum > 80%) were represented in four simplified Bayesian networks (Figure 9).
Finally, the fitting efficiency (e.g., pseudo-coefficient of determination R 2 ), of the ML regression models obtained by the Random Forest algorithm, using as regressors the ecological envelopes variables in comparison to the variables selected by the Bayesian network analysis were very similar. For maritime pine presence (Pb95) models, the fitting efficiency was 69% and 66%, respectively. For maritime pine productivity (SI 50 ) models, the fitting efficiency was 14% and 17%, respectively. For eucalypts' presence (Ec95) models, the fitting efficiency was 72% and 63%, respectively. For eucalypts' productivity (SI 10 ) models, the fitting efficiency was 34% in both cases. These results prove how robust the ecological envelopes variables are in defining species potential areas (Table 6).

Discussion
The results obtained in this study revealed that: (1) the species ecological envelopes indicate that the country has presently a high potential for both species afforestation (over 95%); future climate scenarios forecasted for both species a loss of country's area potentiality, intensified from best to worst scenario (RCP4.5 and RCP8.5), respectively: maritime pine (98%) decreasing to 63-59% in 2050 and to 54-51% in 2070 while eucalypts (95%) decreasing to 73-68% in 2050 and to 62-59% in 2070; (2) it was verified that 99% of maritime pine forest and eucalypts forest in 1995 were found inside excellent, good and regular ecological envelope classes, noting that 62% of maritime pine forest is inside the excellent class compared to the 52% of eucalypts forest; (3) the productivity maps showed that the highest productivity areas were found overlapping areas from regular to excellent ecological envelope classes, namely, 68% for maritime pine and 56% for eucalypts; and (4) the Bayesian network analysis highlighted the most effective variables explaining better the two-species presence in 1995 and their productivities. However, when comparing these selected variables and the ecological envelopes variables as regressors in ML regression tree models the fitting efficiency was very similar proving how robust the ecological envelopes methodology is. The variable thresholds used in species ecological envelopes mapping were consistent with the ones obtained in species distributions (COS1995) and species productivity (1990-1994 plots).
So, despite Portugal having currently a high potential for maritime pine and eucalypts afforestation (over 95%), it is expected that within the next 45-65 years, excellent areas for these species will be limited to a coastal strip in the north of the country, though wider for eucalypts. Moreover, it is expected that present excellent areas will change its potentiality, becoming in the future only good and regular areas for the species. Hence, marginal areas for the species will increase mainly in the inland of the country: maritime pine, from 0.2% to 19-25% and 31-36%, and eucalypts, from 1% to 9-11% and 14-16%, respectively in both scenarios (RCP4.5 and RCP8.5).
Additionally, it is worthwhile noting, when comparing the species present distribution in 1995 with the species ecological envelopes classes, that maritime pine forest was inside the excellent class, while eucalypts new afforestation has been made inside marginal classes for this species. This trend had already been observed in the central inland region of Portugal [67]. This study has also foreseen that this trend will be aggravated in future climate scenarios. Therefore, recommendations on species afforestation potential areas will be key to set back this trend. Additionally, it is important to study the impact of climate change in other Mediterranean species to be used alternatively in future afforestation.
Other studies have reported the decrease of maritime pine and eucalypts potential distribution areas in climate change scenarios, particularly in the south of the Tagus river and in the inland [2][3][4]18]. Overall, it suggested a future trend of species migration from south to north and from the interior to the coastal areas. The eucalypts and maritime pine forests productivity in the northern region may increase in future global change scenarios, mainly in sites closer to the sea and with good quality soils. In the central region, a decrease in those species' productivity will be expected, while in the southern region it is where these two species are expected to be most affected and possibly becoming residual species [2][3][4]19].
The species ecological envelope maps and the productivity maps produced in this study are two approaches to obtain species productive potential. In this study, those maps were produced considering the entire country and with a higher resolution than the ones from previous studies [4,19,22,41], as more input data were used.
Finally, the machine learning approach proved to be a powerful tool to assess the most influential climatic, topographic, and soil variables explaining the two-species distributions and productivities. The results from this study were consistent with previous knowledge about the species [7,34] and with other authors' studies [17][18][19]. Indeed, it was found that the distribution of maritime pine was mainly determined by precipitation-related variables (Figures 1a and 3d), but the elevation and temperature-related variables were very important to differentiate species productivity (Figures 3a,c and 8a). For eucalypts, it was found that this species distribution was mainly explained by temperature-related variables, as well as the species productivity (Figures 1b and 2a,b).
In the light of these findings, it is now understandable why future maritime pine area is forecasted to have a higher decrease in comparison to the eucalypts area. The decreasing precipitation will constrain maritime pine excellent areas to a coastal strip in the north of the country. On the other hand, the increasing temperature will constrain eucalypts excellent areas also to a coastal strip in the north of the country, although wider than for maritime pine, because this temperature increase will allow the species to spread to regions that will become less cold in the future. Although, this temperature increase will also have an opposite effect by turning drier in some regions and constraining this species' distribution area. Additionally, it is important to study the impact of climate change scenarios in other Mediterranean species that may alternatively be used in future afforestation in these areas.
Finally, further investigation is already being undertaken on fitting statistical ML models for the species distribution (SDM) and productivity (SPM) under climate change scenarios to evaluate the statistical precision and agreement of the maps produced in this study (i.e., the ecological envelopes maps and the productivity maps) to the ones obtained by ML modeling.

Conclusions
In this study, the species spatial distribution and productivity maps were obtained for the two main Portuguese wood-production species (maritime pine and eucalypts) for the entire country and with higher resolution than others from previous studies. The species ecological envelope methodology was found to be robust and allowed mapping species potential areas in the present and future projections under climate change scenarios. The machine learning approach allowed understanding the most influential environmental variables that may explain current species distribution and productivity. The distribution of maritime pine was found to be mainly determined by precipitation-related variables, but the elevation and temperature-related variables were very important to differentiate species productivity. For eucalypts, it was found that this species distribution was mainly explained by temperature-related variables, as well as the species productivity.
Understanding which environmental variables will have a higher impact on this two-species distribution and productivity, are key to support recommendations for future afforestation and will bring information for policy design under climate change scenarios. Yet, ecological systems are complex, and many other variables may influence species distribution and productivity. Therefore, further investigation is needed on modeling current and future suitable habitat and productivity by enlarging both the study area and the set of variables studied (for instance, including SoilGrids250m data) to refine species ecological behavior.

Conflicts of Interest:
The authors declare no conflict of interest. Table A1. Summary statistics: Species presence in COS 1995 (maritime pine-Pb95; and eucalypts-Ec95); species productivity in 1990-1994 plots (maritime pine-SI 50 ; and eucalypts-SI 10 ) and the correspondent climatic , topographic and soil variables.