Potential Elevation Shift of the European Beech Stands (Fagus sylvatica L.) in Serbia

According to climate projection models, the global temperature is expected to rise by at least 1.5°C by the end of this century. According to some studies the expected rise in Serbia is even higher. Global warming may result in creating new areas for forest growth. Although creating new forests would be a positive outcome in some areas, global warming can cause negative impacts in other areas, and this can lead to forest loss and the shift of geographical ranges, or even extinction, of plant species. The European beech is the dominant forest tree species in Serbia, featuring high ecological importance and economic value. In mixed or pure stands, beech forests cover approximately 660,400 ha, accounting for 29.3% of the total Serbian forest area. In the present study, the effects of climate change on the distribution of the European beech stands in Serbia, with an emphasis on their elevation shifts, were examined using species distribution models (SDMs). Data for the present tree cover in Serbia, climate projections, and environmental data were used for model building. The models were first tested against present inventory data. In these tests, the models were found to provide accurate projections, as shown by their true skills statistics (TSS) values ranging from 0.652 to 0.736 and area under the curve (AUC) values ranging from 0.868 to 0.937. The potential distribution patterns predicted by the models indicate that the European beech elevational distribution in Serbia would decrease, exhibiting a significant upward shift in elevation during the first part of this century. Current beech stand locations could be changed, and other areas at higher elevations may be more suitable for beech growth. After 2071, European beech stands at elevations below 500 m would be even smaller. This change is caused by temperature rise and occurrence of climate extremes. However, on the highest elevations, further upward shift of the species is not expected.


INTRODUCTION
Over the last 50 years, climate change has been affecting forest ecosystems globally, and climate projection models argue that its impact will increase by the end of this century (Vose et al., 2012;Grimm et al., 2013;Brandt et al., 2016). The expected rise of the global temperature will create new forest areas, particularly in northern zones and at higher elevations. Moreover, increased carbon dioxide concentrations in the atmosphere will further accelerate forest growth (Silva et al., 2016). Although creating new forests is a beneficial outcome of climate change in some areas, other areas would, in contrast, be exposed to high temperature extremes, drought, wildfire, etc., leading to forest loss and the shift of geographical areas or even extinction of plant species (IPCC and Core Writing Team, 2014;Tian et al., 2016). In order to preserve forest ecosystems, many countries are developing strategies and policies to reduce the risks caused by climate change and to seize the opportunities arising from climate change adaptations (Bosworth et al., 2008;Joyce et al., 2009;Littell et al., 2011;Janowiak et al., 2014). Such adaptations are intrinsically related to the management of forest ecosystems, encompassing the determination of climate change impact on an area; to the estimated species' or ecosystem's sensitivity to the projected impacts; as well as to adaptation strategies and their incorporation in the forest ecosystem (Cross et al., 2012(Cross et al., , 2013Stein et al., 2014;Brandt et al., 2016).
Species distribution models (SDMs) are useful tools for developing strategies and adaptation policies to climate change. SDMs create a statistical relationship between the current climate and the occurrence of the species. They use a number of statistical approaches to find relationships and gridded environmental data to extrapolate model projections in space and time (Elith and Leathwick, 2009;Franklin and Miller, 2010;Thuiller and Münkemüller, 2010;Serra-Diaz et al., 2012;Hijmans and Elith, 2017).
The European beech (Fagus sylvatica L.) is the dominant deciduous tree species in Central Europe and in the higher elevated areas of Southern Europe (Bolte et al., 2016;Mausolf et al., 2018). Owing to its high ecological importance and economic value, the European beech has a significant role in the European forestry sector and the overall European forest biodiversity. According to Banković et al. (2009), beech forests, in both mixed and pure stands, cover approximately 660,400 ha in Serbia, accounting for 29.3% of the total Serbian forest area. As the dominant deciduous tree species, it is widely distributed in Serbia. It can be found under different climate conditions (mostly mild winters and moist summers) and at a wide variety of sites with the exception of extremely dry soils featuring a low water storage capacity, stagnic soils, or soils prone to flooding and high groundwater table (Ellenberg and Leuschner, 2010;Bolte et al., 2016). In the future, the increasing variability of climate and the frequency of climatic extremes (IPCC and Core Writing Team, 2014) will impact tree growth (Easterling et al., 2000;Anderegg et al., 2015). A number of studies (Köcher et al., 2009;Zang et al., 2014;Zimmermann et al., 2015;Kunz et al., 2018) consider the European beech as a species sensitive to climatic extremes, especially drought and water deficit, which reduces its competitive advantage over less drought-sensitive species, and this will ultimately result in landscape transformation (Scharnweber et al., 2011;Barigah et al., 2013;Urli et al., 2013). Consequently, the geographical range of the European beech is expected to decrease in the future (Thurm et al., 2018), accompanied by an upward elevation shift of the species. At the present, higher elevations (exceeding 1,200 m) are not favorable to beech growth due to lower temperatures and late frosts (Budeanu et al., 2016;Kolář et al., 2017). According to the national reports on the climate and climate change in the Republic of Serbia (DRINKADRIA project, 2014), the average temperature in Serbia is expected to rise by approximately 3.7 • C by the end of the century (the A2 climate change scenario), so it is expected that beech growth will be viable at higher elevations. However, certain species, such as Norway spruce, will be adversely affected by the forecasted climate change and will become sparser at such elevations (Falk and Hempelmann, 2013).
In the present study, we hypothesize that the European beech area shifts upward in elevation under changing climate conditions. SDMs were used to test the hypothesis in Southeastern Europe in Serbia.

Study Area
The area of this study was the territory of Republic of Serbia situated in Southeastern Europe, in the center of the Balkan Peninsula and the southern part of the Pannonian Plain (between 41 • 53 and 46 • 11 latitude north and 18 • 49 and 23 • 00 longitude east). Available data in the National Forest Inventory used in the study covers 77,474 km 2 , while forests cover amount 22370.17 km 2 , i.e., 25.3% of the territory of Serbia. 1

Presence/Absence Data
Records of the presence/absence of the European beech (F. sylvatica L.) in Serbia were obtained from the National Forest Inventory of Serbia (Figure 1; Banković et al., 2009). These records contain the data for 19,371 tested cells distributed in the 2 × 2-km network (grid), out of which 5,852 were identified as forest cells. Beech was found in 1,651 forest cells, i.e., in an area covering more than 28% of the Serbian forest area. Therefore, beech can be considered as the dominant forest species in Serbia. Beech forests are found in the mountainous regions of Serbia in both pure stands and mixed stands with broadleaves and conifers. Their distribution ranges from 100 to 1,700 m.
The data set was split randomly into 80% training and 20% evaluating subsets. The splitting was stratified in order to keep the data on the presence/absence of the European beech in Serbia equal in the training and the evaluating data subsets. The splitting was carried out using the caret R package (Kuhn, 2008).

Environmental Data Sets
A set of 29 bioclimatic variables was used for both the current and the projected future climate ( Table 1). Bioclim variables were downloaded from the CliMond Archive (Kriticos et al., 2012). The bioclimatic variables were based on the monthly temperature and rainfall values. Those values were used in order to generate biologically meaningful variables, which are relevant for the distribution of beech as they represent annual trends (e.g., the mean annual temperature), seasonality (e.g., the annual range in temperature and precipitation), and extreme or limiting environmental factors. A resolution of 10 arcmin, which FIGURE 1 | Distribution of European beech forests in Serbia according to the National Forest Inventory data (Banković et al., 2009) obtained in R packages "RgoogleMaps" (Loecher and Ropkins, 2015) and "mapproj" (McIlroy, 2018). corresponds to 18.6 × 18.6 km, was employed for the variables, as well as the MIROC-H A2 climate projections.
In order to develop high-precision distribution models, it is important to include soil characteristics in the modeling (Seynave et al., 2008). In addition to the climatic variables, a set of 38 soil variables were used in the study ( Table 1). The first two soil variables describe the soil depth, and for the other soil characteristics, four separate variables were used for the corresponding four soil depths ( Table 1). All of the soil variables used were obtained from the soilgrids.org database using a resolution of 250 m (Hengl et al., 2017).
Elevation maps were derived from the Copernicus Land Monitoring System at a resolution of 25 m. The catdes function of the FactorMineR package (Lê et al., 2008) was used for characterizing variable importance for determining the presence/absence of beech. FactoMineR is R package for multivariate exploratory data analysis. It performs classical principal component analysis (PCA), correspondence analysis (CA), multiple correspondence analysis (MCA), clustering, as well as advanced analyses including different data structures. The model to be used in assessing the climatic change impacts and soil characteristics on the presence of the beech was designed with variables in which p was less than 0.01 (p < 0.01) (Husson et al., 2010). Both single and ensemble SDMs were employed in the study. The ensemble platform biomod2 ) offers the possibility of running 10 modeling techniques for the species distribution modeling used to predict the potential distribution of the European beech. To enhance computation efficiency and minimize computation time requirements, the following seven algorithms were applied: generalized linear models (GLMs), generalized boosting model (GBM), classification tree analysis (CTA), artificial neural networks (ANNs), flexible discriminant analysis (FDA), multivariate adaptive regression spline (MARS), and random forest for classification and regression (RF). Although certain models such as the generalized additive model (GAM), MAXENT.Phillips, and Maximum Entropy MAXENT.Tsuruoka produce good results, they were disregarded in our study for being overly time-consuming in our case and featuring high computer performance requirements.
Owing to substantial discrepancies between the results obtained with different models, the ensemble model was  (Hutchinson et al., 2009;Kriticos et al., 2014;Hengl et al., 2017).

Bio04
Weekly temperature seasonality (C of V)

Bio05
Maximum temperature of the warmest week ( • C)

Bio06
Minimum temperature of the coldest week ( • C)

Bio08
Mean temperature of the wettest quarter ( • C)

Bio09
Mean temperature of the driest quarter ( • C)

Bio10
Mean temperature of the warmest quarter ( • C)

Bio11
Mean temperature of the coldest quarter ( • C)

Bio13
Precipitation of the wettest week (mm)

Bio14
Precipitation of the driest week (mm)

Bio15
Weekly precipitation seasonality (C of V)

Bio16
Precipitation of the wettest quarter (mm)

Bio17
Precipitation of the driest quarter (mm)

Bio18
Precipitation of the warmest quarter (mm)

Bio19
Precipitation of the coldest quarter (mm) considered as the most suitable for minimizing the limitations of single models and producing high-precision results (Coetzee et al., 2009). The performance of each model was assessed using receiver operating characteristic (ROC) known as area under the curve (AUC) and true skills statistic (TSS) values. TSS values were employed because in it, at the range from -1 to 1, either omission or commission errors are not affected by prevalence of the finding under consideration, as opposed to Kappa values (Allouche et al., 2006). According to Gama et al. (2016), AUC values can be random (<0.5), poor (0.5-0.7), good (fair) (0.7-0.9), and excellent (>0.9). Gama et al. (2016) also labels TSS values as poor (0.2-0.5), useful (0.6-0.8), and good to excellent (>0.8).

Test of the Models
For the present climatic conditions, the models showed a good fit between the observed and the predicted presence/absence ( Table 2). The TSS values obtained ranged from 0.652 to 0.736, which is considered as good. A better model performance was indicated by the AUC values ranging from 0.868 to 0.937, which is considered as excellent. Among the single models, RF, GLM, and GBM exhibited the best performance ( Table 2). However, the fundamental feature of the biomod2 package is the ability to combine the projection of single models and build ensemble models based on the single models. The results obtained show that the AUC of the present ensemble model (ENS) was 0.923, which is a bit lower than that obtained using the RF model, but still better than the results obtained using the other models ( Table 2).

Projections for the Periods 2041-2070 and 2071-2100 Under the Climatic Scenario A2
The projected distribution of the European beech stands in Serbia, obtained using the ensemble model and the biomod2 package, indicated that the current beech forest cover will change in the future. Elevational distribution of the beech will be narrowed, and it will be shifting from lower to higher elevations. During 2041-2070, in lower elevations, an upward elevation shift is expected, accompanied by the increased density of the European beech stands at higher elevations (Figure 2). It can be noticed that the maximum density of the beech will be significantly higher in the future. Elevations between 800 and 1,100 m will be more suitable for beech growth than the elevations at which the European beech stands currently predominate in Serbia, i.e., 500-900 m (Figure 2).
After 2071, the projections obtained show that the European beech stands in Serbia found at lower elevations will become even The models were tested with their projections for the presence/absence of European beech in Serbia. For the acronyms of the models, see section "Materials and Methods." Frontiers in Plant Science | www.frontiersin.org sparser, especially at elevations lower than 500 m (Figure 2). The occurrence of beech stands at higher elevations is predicted to be similar as in the previous period (2041-2070), and the species will not shift further upward in elevation. At the highest elevations exceeding 1,400 m, it seems unlikely that beech will migrate above the current upper distribution limit (Figure 2).

DISCUSSION
The results of the present study show that the biomod2 platform is able to predict the current distribution of beech in Serbia. This suggests that it can be used also for predicting any future distribution changes caused by climatic change. The biomod2 platform has been used for such projections in a number of studies (Duque-Lazo et al., 2016;Takolander et al., 2019). In our study, this platform exhibited higher accuracy with regard to single models when tested against the beech distribution in the current climate. The AUC accuracy ranged from 0.868 to 0.937, considered as excellent, whereas the TSS accuracy ranged from 0.652 to 0.736, considered as good. Despite the high accuracy of the single models, the ensemble model forecasting is an effective method for aggregating single SDMs and resolving the problem of intermodel variability, thus leading to results of higher precision (Xu et al., 2015;Lei et al., 2017). The model projections of the present study indicate that European beech areal in Serbia will probably change in the future, mostly in elevation, which will result in landscape transformation of mountain regions. Not only will temperature rise in the future in Serbia, but also the average precipitation is expected to decrease by 15%. The number of tropical and dry days (Kržič et al., 2011) will also increase, which may lead to the shift or extinction of species (Tian et al., 2016). Owing to climate change, many tree species are predicted to shift their geographic ranges (Takolander et al., 2019). As beech is very sensitive to drought, temperature extremes, and global warming, a number of studies (Scharnweber et al., 2011;Kunz et al., 2018) have suggested that extreme climate events, especially drought, can lead to beech forest loss. The results obtained in this study predict an upward beech elevation shift in lower altitudes. The beech shift to higher elevations and beech forest loss have been predicted in many studies (Dorado-Liñán et al., 2017;Sedmáková et al., 2019). Peñuelas et al. (2007) argued that the European beech will shift from lower to higher altitudes in some areas in Spain, accompanied by a decrease in its distribution range. Furthermore, it will replace heathlands and grasslands, as well as some coniferous forests, at higher elevations (Peñuelas et al., 2007). Kolář et al. (2017) concluded that even though average climate characteristics in the montane belt will be, in the future, more favorable for beech, its growth and vitality will be probably still limited by climate extremes, such as late frosts.
In spite of the shift in elevation, it seems unlikely that European beech in Serbia will migrate above current upper distribution limit at elevations exceeding 1,400 m during 2071-2100. This limitation of shifting can be attributed to the sensitivity of beech to drought and temperature rise (Kolář et al., 2017) as well as to shallow mountain soils with a low water storage capacity, which are usually present at higher elevations (Reif et al., 2017). The results obtained in the present study are consistent with the results of Seynave et al. (2008) and Calvaruso et al. (2017), who reported that beech growth is reduced by unfavorable soil properties, especially in shallow soils. However, an upward shift in elevation of the European beech is to be expected. The reason for such an elevation shift is probably the annual temperature rise. Mátyás et al. (2010) reported that a temperature rise of 1 • C can lead to an upward species shift of approximately 170 m along the mountain slope. By comparing pictures taken at three time instants (1 year in the 1920s, 1943, and 2003), Peñuelas et al. (2007) demonstrated that both shifting of beech stands and changes in their density took place simultaneously with warming during the period examined. Yao and Zhang (2015) compared two mountain localities on the Tibetan Plateau with a difference in temperature and showed that the treeline occurred 500-1,000 m higher in the localities with a higher temperature. Moreover, with a temperature rise at higher elevations, soil decomposition will be facilitated (Gutiérrez-Girón et al., 2015) and more nutrients will be available for plant uptake. All these phenomena, accompanied by the phenotypic plasticity of the European beech, can lead to the shift of the species from lower to higher elevations.

CONCLUSION
Using the biomod2 platform, we evaluated the performance of SDMs. As main highlights we found the following: -Models showed good accuracy, and results of the TSS and AUC of the models ranged from good to excellent. -In the middle (2041-2071) and late 21st century (2071-2100), it is expected that habitats of the European beech, which are at elevations between 900 and 1,100 m, will be more suitable for growth than the habitats between 500 and 900 m, at which the beech stands currently predominate in Serbia. -A further upward shift of the beech at the highest elevations above 1,400 m is not expected in Serbia.

AUTHOR CONTRIBUTIONS
LP and DS conceived the presented idea. LP developed the theory and performed the computations. DS and SO verified the analytical methods. DS, ML, and EM encouraged LP to investigate and supervised the findings of this work. All authors discussed the results and contributed to the final manuscript.