Effect of climate change on potential distribution of Cedrus libani A. Rich in the twenty-first century: an Ecological Niche Modeling assessment

The present study is focused on the potential distribution of the Lebanese cedar (Cedrus libani) in the present and in the future throughout the twenty-first century. The location of this work encompasses Lebanon, Syria and Turkey. Twenty-four environmental variables are used and two Representative Concentration Pathways (RCP) scenarios for two different time periods are studied: RCP 4.5 2050, RCP 4.5 2070, RCP 8.5 2050 and RCP 8.5 2070. The most interesting novelty is the use of 13 General Circulation Models and 6 algorithms (Climate Space Model, Envelope Score, Environmental Distance, Genetic Algorithm for Rule-set Production, Maximum Entropy and Support Vector Machines) were considered for modelling. Area Under the Curve is used as goodness of fit and building the final consensus map. The global habitat suitability area would enlarge in the forecasted scenarios with respect to the present, although it would be more restricted in 2070 due to the altitudinal shift. This study also suggests an interesting approach to manage C. libani stands by means of afforestation programs aiming to face global warming in the late twenty-first century.

Introduction summer (Senitza 1989;Atalay 2002). Although average annual rainfall in the different cedar stands vary from 540 up to 1500 mm, the growth is highly affected by spring precipitation patterns and snowmelt (Hajar et al. 2010;Ducrey et al. 2008;Touchan et al. 2005). Cedar tolerates a wide range of temperatures, ranging from − 15 to 7.5 °C with extremes of − 35 °C in winter, and 40 °C in summer (Ducrey et al. 2008).
In Lebanon, the cedar can thrive on alkaline and acidic soils, with no striking difference. It may also be found in karstic-dolomitic rock formations (Abi Saleh et al. 1996). Mainly, C. libani prefers carbonate soils with pH between 6.6 and 8.2 (Ayasligil 1997) but also grows on soils with parent materials consisting of sandstones, mica schist, serpentine and olivine basalt (Sevim 1955). Due to its relative wide range of distribution and its relative plasticity towards climate extremes, coupled with its importance as a major species for timber production, Messinger et al. (2015) investigated its introduction into Central Europe as a promising species for timber production and for ecosystem services, under a changing climate.
Cedar has been widely used in reforestation activities in Turkey and Lebanon, as a major species for ecosystem restoration and for the increase of forest cover (Khuri et al. 2000;Boydak 2003;Ayan et al. 2017). Many programs were devoted to extend the occupancy area of cedar outside its natural range, especially in Central Anatolia (Oner and Uysal 2009). Those actions, which required investments in reforestation by planting and sowing with carpels, doubled the cedar plantation area in Turkey. In Lebanon, the national master plan for land use has foreseen a "cedar and fruit trees corridor" on the western slopes of Mount Lebanon, between 1400 and 1900 m.
In order to optimize these initiatives, it is important to understand the bioclimatic niche and distribution of C. libani at present and under future climate. Such projections might be achieved after a full data collection of present species distribution, which is of paramount importance to clarify the conservation status and the degree of threat the species is facing. Applying ENM (Ecological Niche Modeling) can provide further insights into species' potential distribution and offer the possibility for foresters to plan forthcoming projects devoted to extend the actual distribution, identify new areas where cedar would be able to grow and identify current areas where the species might suffer from changes in climate conditions (Guisan and Zimmermann 2000;Hernández et al. 2006;Pearson et al. 2007;Booth et al. 2014). Correlative approaches are the base for developing ENM, where field observations and environmental variables are used to produce statistically-derived response surfaces. Outputs represent a model that approximates the species' distribution, and/or mainly based on analysis of its natural distribution, but also including some introduced stands outside its natural distribution (c.f. the realized and fundamental niches sensu Hutchinson 1957). This work has been carried out considering 13 GCMs and 6 algorithms.
In view of this, this paper aims to: (1) provide an approximate but reliable distribution of actual cedar stands in the Near East considering political instability, (2) project and map its future habitat suitability using ENM based on environmental variables, and (3) evaluate the spatial distribution of reputed best areas to increase cedar's distribution by new plantations.

Study area
The study area is located in the Near East including Lebanon, Syria and Turkey (Fig. 1a). All the existing populations are within the study area, mainly on the coastal ranges of Lebanon and Syria, between 1300 and 1950 m. Most of the populations are on the western slopes, facing the sea, with few on the leeward established direction. Apart from this natural distribution area, there are many artificial cedar stands established mainly through planting in Aegean, Marmara, Black Sea and Inner and Eastern Anatolian regions of Turkey. They were included in this study because they thrive and show good health. Orography is rough within the area of distribution of the cedar. Mount Taurus range culminates at 3916 m a.s.l. in Erciyes peak, Mount Lebanon reaches 3088 m a.s.l., while Nusayriyah mountain in Syria barely summits at 1575 m a.s.l.

Input data
Twenty-four environmental variables were considered for C. libani occurrence (Table 1). Among them, 3 explanatory variables were topographic (elevation, aspect and slope). A Digital Elevation Model (DEM) at 70 m resolution, from ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) GDEM available at https ://aster web.jpl. nasa.gov/gdem.asp, was used to compute both aspect and slope (30 arc-seconds) by the 3D Analyst function implemented in ArcMap 9.3.1. We used WGS84_UTM_Zone 34 N projected coordinate system. The remaining 21 variables were: (1) 19 BIOCLIM climatic factors downloaded from WorldClim 1.4 project (Hijmans et al. 2005) plus (2) 2 indexes (Emberger's aridity index and Continentality) which were used in the modeling process because these indices were not included in the WorldClim's dataset. They were calculated from the relative formulas as follows: where P is the annual rainfall, M 2 is the mean of the maximum temperatures of the warmest month, and m 2 is the mean of the minimum temperatures of the coldest month (López-Tirado et al. 2018). Conrad Index of Continentality (Ci) was retrieved from the difference between the mean temperatures of warmest (M) and coldest months (m), divided by the sine of the latitude (φ) of the centre of each cell forming the study area (Conrad 1946). The formula is: Raw data in TIFF (Tagged Image File Format) were converted into ASCII (American Standard Code for Information Interchange) format to be used in the modeling software openModeller Desktop v1.5.0.
Occurrence data of C. libani in the Near East were assembled from different databases, personal observations and interpretation of aerial photography. Specifically, populations in Lebanon were retrieved from the forest map MOA (2005), checked and updated by field trips and photo-interpretation. We inspected all polygons drawn into the forest map to validate the attributed land cover and forest type. Additional information, such as centroids' coordinates, polygon surface in hectares and site toponyms, was available in the database. Turkish distribution of cedar was also collected from OGM (2015) as polygons with metadata about surfaces, coordinates and toponyms. Cedar forests in Syria were the hardest to collect and verify, due to known recent political turmoil in that area. Previous investigations made the starting point of an accurate aerial photo-interpretation by satellite images. Two types of images were used: Landsat 5 and IRS-1D (about 5 m of resolution) and images from the General Command of Mapping taken in 2010-2012 (45 cm of resolution). Given the impossibility to verify the real extent of each cedar forest in Syria, such an a priori reconstruction might be considered imperfect, thus over or underestimating the true species occurrence in that country. However, we believe the number of polygons collected in this investigation significantly approximates the distribution in Syria. We did our best to approach the distribution of C. libani in the study area.
Polygons were gathered and stored into a geodatabase to work with ArcMap 9.3.1. On the other hand, a shapefile of polygons was created from a random explanatory variable raster at 30 arc-seconds resolution. Then C. libani polygons were intersected with the vector grid by means of select by location tool of this software. The selected grids were saved as a new polygon Shapefile, which was the base to build a point (centroid) Shapefile using shapes to centroids from XTool Pro extension. A total of 16,363 equidistant points were retrieved as the current distribution of C. libani in the study area. Finally, coordinates associated to each point were then saved in a CSV (Comma-Separated Values) file to be submitted to the modeling software. (EnvScore), Environmental Distance (EnvDist), Genetic Algorithm for Rule-set Production (GARP), Maximum Entropy (MaxEnt) and Support Vector Machines (SVM). Other algorithms were considered at the beginning such as Niche Mosaic, Bioclim and Environmental Niche Factor Analysis (ENFA), but then dismissed because of non-reliable results in AUC values. A consensus map across the six algorithms was built following the equation (Vessella et al. 2015(Vessella et al. , 2017López-Tirado et al. 2018):

Ecological Niche Modeling and goodness of fit
where X pond is the weighted probability of habitat suitability, z i is the probability of the potential occurrence of the cell of the i-model and AUC i is the Area Under the Curve. AUC values were also used to validate the models according to Swets (1988) who established the following intervals: 0.5-0.7 as low accuracy; 0.7-0.9 as potentially useful; and > 0.9 as high accuracy.
The final consensus map for the present plus the expected future explanatory variables downloaded from WorldClim 1.4 project (Hijmans et al. 2005) were the base for building the scenario maps along the twenty-first century-reclassified in ten intervals by means of ArcMap 9.3.1. The latter were performed in two RCPs (Representative Concentration Pathways), each one in two time slices. These scenarios represent the concentration of greenhouse gases and pollutants resulting from human activity, adopted by the Fifth Assessment Report (IPCC 2014). The most extreme and an intermediate scenario were selected: the former achieving 8.5 W/m 2 by 2100 (RCP 8.5), while the latter an increase of 4.5 W/m 2 (RCP 4.5), also in the same deadline. Two time slices were studied for each RCP:

Results
The total study area spans ca. 978,202 km 2 within which the observed occurrence of C. libani only reaches 0.39%. Our results suggest a greater potential distribution of this species in the present reaching the 9.46% of the whole territory in the consensus map for those cells above 50% of suitability. This situation could be explained by habitat fragmentation. The habitat suitability area for the future scenarios is increased as well ranging from 21.04 in RCP 8.5 (2070) to 41.66% in RCP 4.5 (2050) (see Table 2 for the results in km 2 and %). The potentially suitable habitat of C. libani in the present occurs along Southern Turkey, where most of the stands are found today, and not reaching the top of the mountains (Fig. 1b). Figure 2 shows the output maps for the RCPs scenarios. All the scenarios follow similar trends with respect to Fig. 1b: (1) a greater habitat suitability area, and (2) a higher elevation and latitude (cf. Table 2). As can be noted in Fig. 1b the maximum suitability value reached was 80%. In other words, no cells with 100% for each algorithm were found. Table 2 Elevation, latitude and habitat suitability area for the current distribution and models in the present and forecasted scenarios Values for the potential distribution in the present and forecasted scenarios were computed with cells whose suitability is either equal or over 50% According to our results, Central-Northern Turkish ranges become more suitable for the growth of C. libani in the forecasted scenarios. On the other hand, Southern Turkish ranges lose potential distribution at lower elevations whilst summits are reached by this conifer. The total habitat suitability area is decreased gradually from the intermediate emission scenario (RCP 4.5) in 2050 to the high emission scenario (RCP 8.5) in 2070 (Table 2). Regarding Syria and Lebanon, the main potential stands of C. libani are found in coastal ranges that are influenced by the Mediterranean climate. Inner areas of Syria do not exhibit any potentiality due to the arid conditions-poor occurrence data must be taken into account in this country.
The suitability variation map among the 6 algorithms used to build the consensus map for the current time period is shown as supplementary material ( Figure S1). Weighted Standard Deviation (S.D.) has been used scoring maximum values between 0.41 and 0.44% in some coastal areas of the whole territory and inner areas of Turkey. Goodness of fit of the 6 algorithms was assessed by means of the AUC scores. EnvDis and SVM were the most accurate algorithms whereas EnvScore and MaxEnt showed lower values (Table 3). Figure 3 shows the surface occupancy in the present and both emission scenarios and periods respectively. The four studied scenarios follow a similar tendency where peaks (ca. 20-27% of surface occupancy) are reached at 40-50% of suitability. Between RCP 4.5 and RCP 8.5, only the latter exceeds the 25% of surface occupancy in cells of ca. 40% of suitability. Present period differs from the future scenarios showing lower surface occupancy at this point of suitability. The highest values are found around 32% of suitability. Finally, the highest suitability (over 60%) presents lower surface occupancy (around 0-10%) for the five models.

Discussion
Ecological Niche Modeling has been widely used in the recent decades to understand spatial and temporal potential distribution of many plant and animal species. Most of these theoretical approaches focus on forecasting due to the expected global climate change in the twenty-first century. The Mediterranean Basin would suffer a generalized rising of the mean temperatures whereas rainfall would undergo irregular patterns (IPCC 2014). In addition, Regato and Salman (2008) state that Mediterranean mountainous areas are at the forefront of global climate change (increasing drought effect, rising and fluctuating temperatures, irregular and heavy precipitation, etc.) in the world. Recent studies have demonstrated that mountain species such as Cedrus atlantica are persisting in restricted and isolated areas considered as microrefugia (Cheddadi et al. 2017). Increasing drought effects are harmful to sensitive species like mountain conifers and Cedrus in particular (Linares et al. 2013). Ulbrich et al. (2006) examines various scenarios; in the Mediterranean Basin, towards the end of the century, the net water shortage may increase by 1.3-1.4%, while the annual rainfall may decrease by 10-30%. At this point, we have centered on an emblematic cedar (C. libani) from the Near East which could be affected by the changing conditions in the next decades.
The use of planted stands was justified here because they grow and thrive; therefore they provide information about the species climatic adaptability. Some species, especially trees, could adapt to broader climatic conditions than known in their natural distribution. Thirteen GCMs calibrated for the Northern Hemisphere and 6 different algorithms were used. This methodology has already been proven and achieved previously (López-Tirado et al. 2018). Goodness of fit of the models was executed by the AUC showing disparate scores and only two algorithms displayed a value under 0.7 (Table 3), i.e. indicating low accuracy according to Swets (1988). Nonetheless, they were considered in the consensus map because: (1) their value was close to the threshold (0.68 and 0.64 for EnvScore and MaxEnt respectively) and (2) their weight in the final consensus map was lower according to the equation showed in Materials and Methods section.
Surface occupancy follows similar trends by period (2050 and 2070) in the two emission scenarios. It can be noted from Fig. 3 how present shows the highest surface occupancy at lower suitability class with respect to both future periods. According to these results, C. libani could find wider suitable areas than in the current period. The increase in the area covered by middle suitability (40%) when compared to the present situation could be related to the shifting in latitude in the area of extent of C. libani. The cedar would reach larger areas in the northern part of Turkey under Euxinian climate (a type of climate typical from northern aspects of the Black Sea Region) characterized by a totally  (2050 and 2070); present is also shown different precipitation pattern (substantial summer rain) and higher humidity. As a result, the cedar would be in competition with species that are in their optimal zone (Fagus orientalis Lipsky, Castanea sativa Mill., Carpinus spp., Quercus spp., Abies nordmanniana (Steven) Spach, Pinus nigra J.F.Arnold, etc.) and in suitable conditions for fungus outbreaks. It must be emphasized that this work is addressing the potential distribution of C. libani. Therefore, habitat suitability according to the explanatory variables is detected; no biotic interactions are considered because of their difficulty to be incorporated into ENMs. The suitability of the Black Sea Region and the Central Anatolia regions for the species indicated in the results of the study is confirmed by the performances of the afforestation in the last 30 years in Turkey (Fig. 1a). Whereas, on the southern part of its area of distribution, pest outbreak (i.e. Cephalcia tannourinensis) would be more recurrent and damaging since the insect proliferates under reduced soil moisture and increased soil temperature (Nemer et al. 2014).
Study results indicate that annual rainfall is the best predictor of establishing successful plantations. Where the average annual rainfall falls below 400-500 mm, even if the soil conditions are very favorable, the growth of the cedar plantation is adversely affected (Ürgenç 1986;Karataş and Özkan 2017). This emphasizes that climate characteristics are more effective than the other site characteristics for the height growth of C. libani. In addition, it was determined that there is a positive relationship between height growth and average high temperature, potential evapotranspiration, water excess, average temperature of coldest month, average temperature of warmest month and average temperature of three summer months.

Conclusions
Migration might be necessary for C. libani to cope with projected change. A similar trend has been detected in Cedrus atlantica populations in Algeria where its southern distribution limit shows a high level of tree mortality during recent decades (Slimani et al. 2014). Factors that are correlated with altitude seem to be the most significant (Körner 2007). This is consistent with other tree species studied in the Mediterranean basin (López-Tirado et al. 2018). Habitat suitability area results wider in the forecasted scenarios than in the present period. These results can be useful to manage C. libani natural and non-natural stands. Afforestation programs by ENM models could be carried out to deal with the expected climate change together with other disciplines such as remote sensing (Camarretta et al. 2020).