MaxEnt Modeling for Predicting the Current and Future Potential Geographical Distribution of Quercus libani Olivier

: The purpose of the study was to model the current and potential future distribution of Quercus libani Olivier (Lebanon Oak), a tree species in Turkey, and to predict the changes in its geographical distribution under di ﬀ erent climate change scenarios. In this study, 19 bioclimatic variables at a spatial resolution of 30 arc seconds (~1 km 2 ) were used, collected from the WorldClim database. The bioclimatic data with high correlation according to 31 sets of presence data on the species were reduced with principal component analysis (PCA), and the current and potential distribution were identiﬁed using MaxEnt 3.4.1 software. In order to predict how the distribution of the species will be a ﬀ ected by climate change, its potential geographical distribution by 2050 and 2070 was modeled under the Representative Concentration Pathways (RCP) RCP 4.5 and RCP 8.5 scenarios of the species using the Community Climate System Model (CCSM, version 4), which is a climate change model based on the report of the ﬁfth Intergovernmental Panel on Climate Change (IPCC). Change analysis was performed to determine the spatial di ﬀ erences between its current and future distribution areas. The study results showed that the suitable areas for the current distribution of Quercus libani Olivier cover 72,819 km 2 . Depending on the CCSM4 climate model, the suitable area will decline to 67,580 km 2 by 2070, according to the RCP 4.5 scenario, or 63,390 km 2 in the RCP 8.5 scenario. This may lead to a reduction in the future population of this species. The change analysis showed that suitable and highly suitable areas will decrease under global climate change scenarios (RCP 4.5 and RCP 8.5) for both current and future potential distribution areas. In this context, our study results indicate that for the management of this species, protective environmental measures should be taken, and climate change models need to be considered in land use and forest management planning.


Introduction
Rapid climate change, habitat fragmentation, invasion by foreign species, water, soil, and air pollution, exploitation of nature, and human population growth are the most important factors that destroy the structural and functional integrity of ecosystems [1,2]. Climate parameters (temperature, precipitation, humidity, etc.) have a significant impact on the growth and development of plants, and thus determine the geographical distribution of plant species [3,4]. Climate systems tend to change naturally throughout world history (around 4.5 billion years) [5]. In reaction to these changes, extreme environmental conditions were created that caused the extinction of some species or changed their habitats or niches [6]. People are now aware of global climate change. Plant communities are particularly vulnerable to changes in climate [7,8].
Today, changes in the global climate are unprecedented. Global warming could be more severe than expected and will continue to rise over the 21st century by a minimum of 0.3-1.7 • C to a maximum of 2.6-4.8 • C [9]. In particular, plant species with narrow habitats are under threat due to climate change and are more vulnerable to extinction [10][11][12][13][14][15][16][17]. Therefore, strict precautions should be taken to avoid the permanent loss of endangered species. In this context, suitable habitat prediction and mapping are critical to monitor endangered species, species in danger of extinction, and decreasing native communities [18].
The effect of climate change on endangered plant species in terrestrial ecosystems is complicated, and further information about potential distribution areas is needed to identify and rehabilitate them. The current potential and future distribution, according to different climate change scenarios, can be determined through machine learning methods by means of the layers created using the records of point areas where species are currently, and digital bioclimatic data regarding these areas [19][20][21][22][23][24][25][26][27][28]. Various species distribution models such as CLIMEX and Doup estimate the impact of environmental and climate changes on species and ecosystems. The Genetic Algorithm for Rule Set Production (GARP) produces a rules set, while maximum entropy modeling (MaxEnt) aims to recover critical habitats [29][30][31]. Among these modeling approaches, MaxEnt has advantages because it uses both presence data and categorical data as input variables. In addition, it tests forecast precision, it is always stable and reliable, even in small sample sizes and with gaps, and it directly produces a spatially open habitat suitability map and evaluates the significance level of individual environmental variables using a built-in jackknife test [20,[32][33][34][35].
The main purpose of this study is to predict the current and future potential geographical distribution of Quercus libani Oliv. (Lebanon Oak) in Turkey. Presence data and bioclimatic variables were used and a prediction model was created using a species distribution model with MaxEnt software under the global climate change scenarios RCP 4.5 and RCP 8.5. We explored suitable areas for Quercus libani Oliv. and analyzed the effect of each bioclimatic variable with PCA analysis. Additionally, the changes were mapped by comparing the current geographical distribution of Q. libani Oliv. with its predicted geographical distribution in 2050 and 2070, according to the RCP 4.5 and RCP 8.5 climate scenarios. We aim to produce information about the Lebanon Oak by the study area and to promote a new approach to forestry-related studies. Our research questions are as follows: 1. What is the effect of each independent bioclimatic variable in the creation of the model? 2. What is the current and future potential geographical distribution of Quercus libani Oliv. (Lebanon Oak) in Turkey?
3. What are the suitable areas for Quercus libani Oliv. according to global climate change scenarios for the future (2050 and 2070)?
4. What will the changes be from the current potential geographical distribution of Quercus libani Oliv. according to the RCP 4.5 and RCP 8.5 climate scenarios in 2050 and 2070? Finally, we discuss the possible implications of the results for future global forestry research.

Materials and Methods
Oaks (Quercus L.) are crucial species in global forestry compared to other broad-leaved species, with large forests in the temperate zones of North America, Europe, and East Asia containing more than 200 species, along with several subspecies, varieties, and natural hybrids [36]. The total cover area of oak stand in the forests in Turkey is 5.8 million hectares, accounting for 26% of forest forms, which makes it the predominant tree type in Turkey [37]. Quercus libani Oliv., one of 24 oak species that naturally exist in Turkey, has a local distribution specific to the east Mediterranean, including West Syria, Northeast Israel, East Turkey, North Iraq, and West Iran. It is a semi-evergreen or deciduous forest tree species growing to 8-10 m (Figure 1). It is usually distributed at an elevation of 700-2000 m [38]. With its limited geographical distribution in Turkey, this species is at risk of becoming a threatened or endangered species. It is in the Least Concern (LC) category of the red list of the International Union for Conservation of Nature (IUCN) and has a decreasing population [39,40]. As the wood of the Lebanon oak is resistant to insects and fungi, it is preferred in construction work. Furthermore, it is an important natural material for the pharmaceutical industry due to the high tannin content of its leaves and fruits, and it is used in biological control to keep pests away from plants [41].
Sustainability 2020, 12, x FOR PEER REVIEW  3 of 19   of the International Union for Conservation of Nature (IUCN) and has a decreasing population [39,40]. As the wood of the Lebanon oak is resistant to insects and fungi, it is preferred in construction work. Furthermore, it is an important natural material for the pharmaceutical industry due to the high tannin content of its leaves and fruits, and it is used in biological control to keep pests away from plants [41].  [42] In this study, the coordinates of the points representing the geographical distribution of Quercus libani Oliv. in Turkey selected as the study area were determined and marked on the map ( Figure 2). Thirty-one sample points were identified for the presence data used in the study according to the literature and available data in areas where the species is distributed naturally [43][44][45][46]. The coordinates of these points were marked on a WGS84 coordinate system using high-resolution Google Earth satellite image data and QGIS 3.10.0 [47] software. The descriptive data of the sample points are presented in Table 1.  In this study, the coordinates of the points representing the geographical distribution of Quercus libani Oliv. in Turkey selected as the study area were determined and marked on the map ( Figure 2). Thirty-one sample points were identified for the presence data used in the study according to the literature and available data in areas where the species is distributed naturally [43][44][45][46]. The coordinates of these points were marked on a WGS84 coordinate system using high-resolution Google Earth satellite image data and QGIS 3.10.0 [47] software. The descriptive data of the sample points are presented in Table 1. of the International Union for Conservation of Nature (IUCN) and has a decreasing population [39,40]. As the wood of the Lebanon oak is resistant to insects and fungi, it is preferred in construction work. Furthermore, it is an important natural material for the pharmaceutical industry due to the high tannin content of its leaves and fruits, and it is used in biological control to keep pests away from plants [41].  [42] In this study, the coordinates of the points representing the geographical distribution of Quercus libani Oliv. in Turkey selected as the study area were determined and marked on the map ( Figure 2). Thirty-one sample points were identified for the presence data used in the study according to the literature and available data in areas where the species is distributed naturally [43][44][45][46]. The coordinates of these points were marked on a WGS84 coordinate system using high-resolution Google Earth satellite image data and QGIS 3.10.0 [47] software. The descriptive data of the sample points are presented in Table 1.   The bioclimatic variables in the WorldClim database were used for modeling the current potential distribution areas. WorldClim (version 1.4) contains information about monthly minimum, maximum, and mean temperatures from 1960 to 1990, as well as monthly mean areal precipitation. The bioclimatic variables used to predict the current distribution and having a spatial resolution of 30 s (~1 km 2 ) were derived from the observational data in WorldClim ver. 1.4 [48,49].
The Community Climate System Model (CCSM ver. 4) was used to predict the future geographical distribution of Quercus libani Oliv. in the study. This version simultaneously shows the atmosphere, ocean and land cover, and sea ice as well as the central coupler component in May 2010. It has four models that allow researchers to conduct basic studies about the climate in the past, at present, and in the future. As the climate scenario related to this model, the RCP 4.5 and RCP 8.5 scenarios from the fifth Report of the Intergovernmental Panel on Climate Change (IPCC5) were used. RCPs are used in climate modeling studies and research to describe potential climate scenarios that are thought to depend on greenhouse emissions in the near future [50][51][52]. In future scenarios, the bioclimatic data for 2050 represent the mean values from 2041 to 2060, while the data for 2070 represent the mean values from 2061 to 2080 [53].
Statistical analysis showed that some of the bioclimatic variables of the sample points have a very high correlation. In order to ensure clarity among bioclimatic variables that have shown similarity, a principal component analysis (PCA) was performed. We identified the possibility of the PCA in our study through a Kaiser-Meyer-Olkin (KMO) test. Field [54] stated that the lower limit of the KMO should be 0.5 to factorize a data cluster. The Bartlett test allows for deciding if the data matrix is a unit matrix and the correlations among the variables are sufficient. All correlation coefficients are 0 to test the H 0 hypothesis [55]. KMO was calculated with the primary assessment using 19 bioclimatic variables. The variables that lack clarity were discounted in the final calculation [56]. Following the KMO test, PCA was performed for identifying variables. PCA is used to describe the low number of principal components that effectively summarize most of the variations in the data and reduce the dimensionality of the problem [57].
With the purpose of modeling the potential and future distribution of Quercus libani Oliv., the maximum entropy approach was used with the MaxEnt 3.4.1 software [30,58]. For the model, 25% of the presence data (seven items) were separated out as the test data, and the MaxEnt modeling procedure was applied with the auto features. To determine the performance of the model, values in the area under the curve (AUC) from the receiver operating characteristics (ROC) were analyzed [24,59]. The AUC value obtained can be interpreted as the estimated probability of the randomly selected grid cell in a correctly adjusted model. AUC defines the success of the model with all possible thresholds. If AUC > 0.5, it means that the model performs better than a random estimation [60]. The closer the AUC test value is to 1, the better the separation is, and the more sensitive and descriptive the model is [61]. To interpret the AUC value, the threshold values are as follows: AUC ≥ 0.9 = very good, 0.9 > AUC ≥ 0.8 = good, and AUC < 0.8 = poor [62,63]. Finally, to determine the contribution of environmental variables, the Jackknife test option was applied in the MaxEnt software [32,64]. This test enables us to determine the significance of each independent variable in the creation of the model. Five thresholds were used for the distribution area to create potential distribution maps: 0 means unsuitable, 0-0.25 means not very suitable, 0.25-0.50 means moderately suitable, 0.50-0.75 means suitable, and 0.75-1 means highly suitable.
The geographical distribution maps for 2050 and 2070 created according to two different climate change scenarios were compared with the current potential distribution map. With a view to determining the changes [65], the data from three different time points were converted to polygon data with the raster/vector conversion function using GIS. Intersection analysis was performed for the polygon data to calculate the differences.

Statistical Analysis of Bioclimatic Variables
The variables, apart from the cluster sets created in the data space, were extracted from 19 bioclimatic variables used as the dataset in the study, and the KMO value was calculated as 0.65 for the remaining nine variables ( Table 2). As this value was higher than the limit value required for factorizing the data cluster, a PCA was applied to the data. Regarding the data groups created as a result of the PCA, nine bioclimatic variables forming the first two groups represented 95% of the total variance ( Table 3). The first group in the reliability analysis, which had a Cronbach's alpha value of 0.986 and was comprised of Bio1, Bio5, Bio6, Bio9, Bio10, and Bio11 was named the temperature group. The second group, which had a Cronbach's alpha value of 0.767 and was comprised of Bio12, Bio13, and Bio16 was called the precipitation group. Figure 3 shows the components, rotated according to the varimax method.  The bioclimatic variables selected for the study in the QGIS 3.10.0 were uploaded to MaxEnt software and the current geographical distribution of Quercus libani Oliv. was modeled. In order to test the accuracy of the model, the AUC values of the training data and test data were analyzed (Figure 4). In terms of the AUC values, precise results were obtained for the training data with 0.91 The bioclimatic variables selected for the study in the QGIS 3.10.0 were uploaded to MaxEnt software and the current geographical distribution of Quercus libani Oliv. was modeled. In order to test the accuracy of the model, the AUC values of the training data and test data were analyzed ( Figure 4). In terms of the AUC values, precise results were obtained for the training data with 0.91 and acceptable results were obtained for the test data with 0.83. The model was considered to be sensitive and descriptive. To determine the contribution of the bioclimatic variables, the Jackknife test option in the MaxEnt modeling program were used ( Figure 5). This test revealed that the Bio6, Bio11, and Bio13 variables made a relatively higher contribution in the model. Moreover, all bioclimatic variables used in the model were understood to be necessary. The bioclimatic variables selected for the study in the QGIS 3.10.0 were uploaded to MaxEnt software and the current geographical distribution of Quercus libani Oliv. was modeled. In order to test the accuracy of the model, the AUC values of the training data and test data were analyzed ( Figure 4). In terms of the AUC values, precise results were obtained for the training data with 0.91 and acceptable results were obtained for the test data with 0.83. The model was considered to be sensitive and descriptive. To determine the contribution of the bioclimatic variables, the Jackknife test option in the MaxEnt modeling program were used ( Figure 5). This test revealed that the Bio6, Bio11, and Bio13 variables made a relatively higher contribution in the model. Moreover, all bioclimatic variables used in the model were understood to be necessary.

Current and Future Potential Geographical Distribution of Quercus Libani Olivier and Its Suitability
According to Global Climate Change Scenarios Figure 6 shows the current and future potential distribution of Quercus libani Oliv. under the global climate change scenarios RCP 4.5 and RCP 8.5 for 2050 and 2070. The model was classified into five levels (unsuitable to highly suitable) in the maps. The current and future potential distribution of Quercus libani Oliv. can be tracked on the maps. Table 4 shows the current and future distribution of this species according to the suitability classes on these maps.
According to the RCP 4.5 climate change scenario, there is no notable change by 2050, but there tends to be a decline in the potential geographical distribution of Quercus libani Oliv. by 2070. In the

Current and Future Potential Geographical Distribution of Quercus Libani Olivier and Its Suitability
According to Global Climate Change Scenarios Figure 6 shows the current and future potential distribution of Quercus libani Oliv. under the global climate change scenarios RCP 4.5 and RCP 8.5 for 2050 and 2070. The model was classified into five levels (unsuitable to highly suitable) in the maps. The current and future potential distribution of Quercus libani Oliv. can be tracked on the maps. Table 4 shows the current and future distribution of this species according to the suitability classes on these maps.

Changes from the Current Potential Geographical Distribution of Quercus libani Olivier
The changes were identified by comparing the current geographical distribution of Q. libani Oliv. with its predicted geographical distribution in 2050 and 2070 under the RCP 4.5 and RCP 8.5 climate scenarios (Figures 7 and 8).  According to the RCP 4.5 climate change scenario, there is no notable change by 2050, but there tends to be a decline in the potential geographical distribution of Quercus libani Oliv. by 2070. In the RCP 8.5 scenario, however, it will have a narrowing geographical distribution in 2050 and 2070. These predictions suggest that the geographical distribution of Quercus libani Oliv. will decrease in the future, while the possible habitat losses will be severe according to both scenarios.

Changes from the Current Potential Geographical Distribution of Quercus libani Olivier
The changes were identified by comparing the current geographical distribution of Q. libani Oliv. with its predicted geographical distribution in 2050 and 2070 under the RCP 4.5 and RCP 8.5 climate scenarios (Figures 7 and 8).    With the purpose of determining the changes, the data from two different time points were compared. The change was considered in two dimensions in this study. The spatial change was determined in the first dimension. The geographical distribution of the species was modeled for different time points and the spatial changes were indicated (Figures 7 and 8). A spatial analysis of these changes was planned for the second dimension. To do so, the current potential geographical distribution of Q. libani Oliv. was classified. These classes, named suitability classes, show the potential geographical distribution for the species, from highly suitable ones to unsuitable ones. By using these classes, both the current and future potential geographical distribution according to the RCP 4.5 and RCP 8.5 climate scenarios in 2050 and 2070 were calculated. The changing geographical distribution areas from the present day to 2050 and 2070 were compared on the change maps, the change in the suitability classes was identified, and the surface areas of these changes were calculated. The codes that belong to current and future suitability were therefore compared. Due to the predicted geographical distribution, the direction and size of the changes were calculated (Table 5, Figure 9).   With the purpose of determining the changes, the data from two different time points were compared. The change was considered in two dimensions in this study. The spatial change was determined in the first dimension. The geographical distribution of the species was modeled for different time points and the spatial changes were indicated (Figures 7 and 8). A spatial analysis of these changes was planned for the second dimension. To do so, the current potential geographical distribution of Q. libani Oliv. was classified. These classes, named suitability classes, show the potential geographical distribution for the species, from highly suitable ones to unsuitable ones. By using these classes, both the current and future potential geographical distribution according to the RCP 4.5 and RCP 8.5 climate scenarios in 2050 and 2070 were calculated. The changing geographical distribution areas from the present day to 2050 and 2070 were compared on the change maps, the change in the suitability classes was identified, and the surface areas of these changes were calculated. The codes that belong to current and future suitability were therefore compared. Due to the predicted geographical distribution, the direction and size of the changes were calculated (Table 5, Figure 9).

Discussion
As a result of global climate change, scientists have developed algorithms to provide useful data towards a better understanding the current distribution of plants and to help predict their future distribution areas based on several environmental variables [66,67]. MaxEnt is popular software based on maximum entropy [68][69][70]. In this study, the current and future distribution areas for Q. libani Oliv., a regional species in Turkey, were modeled using MaxEnt. Our results confirmed that the highly suitable area for this species is Eastern Anatolia [37,38,44]. However, in the 2050 and 2070 scenarios, Q. libani Oliv. will disappear if its distribution is limited to the current area of Southern Anatolia.
Bioclimatic variables predicting species distribution areas do not necessarily indicate the importance of topography. Topography plays a crucial role in the distribution of species in terms of its effect on the microclimate [71,72]. Related to the topography, there may be some suitable areas in adjacent countries like Persia, Syria, and Lebanon, but the distribution models for Turkey created in this study cannot provide information on them [73].
Although Q. libani Oliv. is an important oak species, there is limited research on it. Furthermore, we could not find any related research about its geographical distribution area in the recent literature. In this study, the occurrence of Q. libani Oliv. is predicted by determining not only the current geographical distribution, but also the impact of climate change on its future. For this prediction, we used MaxEnt software.
It is clear that the distribution areas and population of Q. libani Oliv. will change, but the nature of that change depends on the climate change scenario. Similar studies conducted on this topic also suggest that climate change is one of the most important factors that will affect habitat losses and fragmentation as well as the loss of biological diversity [16,19,[74][75][76]. In terms of biogeography, such changes and as dislocations, population reductions, and species extinctions are perennial challenges [77][78][79].

Discussion
As a result of global climate change, scientists have developed algorithms to provide useful data towards a better understanding the current distribution of plants and to help predict their future distribution areas based on several environmental variables [66,67]. MaxEnt is popular software based on maximum entropy [68][69][70]. In this study, the current and future distribution areas for Q. libani Oliv., a regional species in Turkey, were modeled using MaxEnt. Our results confirmed that the highly suitable area for this species is Eastern Anatolia [37,38,44]. However, in the 2050 and 2070 scenarios, Q. libani Oliv. will disappear if its distribution is limited to the current area of Southern Anatolia.
Bioclimatic variables predicting species distribution areas do not necessarily indicate the importance of topography. Topography plays a crucial role in the distribution of species in terms of its effect on the microclimate [71,72]. Related to the topography, there may be some suitable areas in adjacent countries like Persia, Syria, and Lebanon, but the distribution models for Turkey created in this study cannot provide information on them [73].
Although Q. libani Oliv. is an important oak species, there is limited research on it. Furthermore, we could not find any related research about its geographical distribution area in the recent literature. In this study, the occurrence of Q. libani Oliv. is predicted by determining not only the current The findings of the study show that the spatial distribution of this species will decrease by around 10% by 2070 according to the RCP 8.5 scenario. The proportion of highly suitable areas is estimated to decrease from 4% to 1.5% in the future.
Almost all bovine and ovine animals feed on oak leaves and acorns. Therefore, areas where oak species are present are considered to be the most suitable places for grazing [80]. The current distribution map demonstrates that this species is distributed in Hakkari,Şırnak, Van, Bitlis, Muş, Bingöl, Malatya, K. Maraş, Tunceli, and Elazıg, as well as Adıyaman, Diyarbakır, Batman, and Mardin provinces where livestock husbandry is the most important agricultural activity. Therefore, the species is under pressure from grazing. Excessive or unplanned grazing will further threaten the existence of this species in the future. For that reason, areas with young individuals and degraded stands should be protected.
There are 341 seed stands in Turkey, only 13 of which belong to oak species [81]. Out of these oak seed stands, there is only one for Q. libani Olivier. This number is insufficient and should be increased for Lebanon Oak. To this end, the habitat conditions of the existing stands should be identified, and the optimum ones retained, with some of them reserved as seed stands.
The RCP 4.5 scenario estimated that there will be no major changes in the geographical distribution of the species by 2050. However, there will be a decline in suitable and highly suitable areas by 2070. The results of the analysis made for the RCP 8.5 scenario show that the losses in terms of its geographical distribution are predicted to be severe. Although the prediction in this scenario does not indicate any marked change for 2050, there will be a major decrease in suitable and highly suitable areas by 2070.
As regards the direction of the changes expected to take place in the geographical distribution of Q. libani Oliv., only 4% of the currently unsuitable areas for its potential distribution will become suitable by 2050 according to the RCP 4.5. scenario, while 26% of the suitable areas will become unsuitable. By 2070, those rates will be 3% and 38%, respectively. By 2050, the rates will be 3.9% and 46%, respectively, and 4.6% and 56%, respectively, by 2070, according to the RCP 8.5 scenario. These rates indicate that the RCP 8.5 scenario envisions more severe consequences.

Conclusions
This study revealed that the geographical distribution areas of Q.libani Oliv. will be negatively affected by climate change according to a model based on the maximum entropy algorithm, and the suitable distribution areas will be narrower. Therefore, we suggest that the state forestry authority, which manages 99% of the forests in Turkey, should take precautions for the future. Moreover, it should be emphasized that land use plans should be developed and genetic studies should be supported to conserve seed sources for the continuity of the species and ensure the species' transition to possible distribution areas in the future, apart from the protective measures put in place to prevent excessive grazing and consumption. Stands with Q. libani Oliv. individuals that are resistant to drought can be identified through genetic studies to ensure the continuity of the species and its resistance to climate change in the future. Therefore, efforts should be undertaken to conserve for the future for those individuals that will be less significantly affected by the negative impacts of climate change and can continue surviving and reproducing.

Conflicts of Interest:
The authors declare no conflict of interest.