Past, present and future distributions of Oriental beech (Fagus orientalis) under climate change projections

RESEARCH ARTICLE Past, present and future distributions of Oriental beech (Fagus orientalis) under climate change projections Dilsad DagtekinID1☯¤*, Evrim A. Şahan1☯, Thomas Denk2‡, Nesibe Ko¨ se3‡, H. Nu¨ zhet Dalfes1‡ 1 Eurasia Institute of Earth Sciences, Istanbul Technical University, Istanbul, Turkey, 2 Department of Paleobiology, Swedish Museum of Natural History, Stockholm, Sweden, 3 Department of Forest Botany, Faculty of Forestry, Istanbul University-Cerrahpasa, Istanbul, Turkey ☯ These authors contributed equally to this work. ¤ Current address: Department of Evolutionary Biology and Environmental Studies, University of Zurich, Zurich, Switzerland ‡ These authors also contributed equally to this work. * dilsad.dagtekin@ieu.uzh.ch Abstract Species distribution models can help predicting range shifts under climate change. The aim of this study is to investigate the late Quaternary distribution of Oriental beech (Fagus orien- talis) and to project future distribution ranges under different climate change scenarios using a combined palaeobotanical, phylogeographic, and modelling approach. Five species distri- bution modelling algorithms under the R-package ‘biomod2‘were applied to occurrence data of Fagus orientalis to predict distributions under present, past (Last Glacial Maximum, 21 ka, Mid-Holocene, 6 ka), and future climatic conditions with different scenarios obtained from MIROC-ESM and CCSM4 global climate models. Distribution models were compared to palaeobotanical and phylogeographic evidence. Pollen data indicate northern Turkey and

the western Caucasus as refugia for Oriental beech during the Last Glacial Maximum. Although pollen records are missing, molecular data point to Last Glacial Maximum refugia in northern Iran. For the mid-Holocene, pollen data support the presence of beech in the study region. Species distribution models predicted present and Last Glacial Maximum dis- tribution of Fagus orientalis moderately well yet underestimated mid-Holocene ranges. Future projections under various climate scenarios indicate northern Iran and the Caucasus region as major refugia for Oriental beech. Combining palaeobotanical, phylogeographic and modelling approaches is useful when making projections about distributions of plants. Palaeobotanical and molecular evidence reject some of the model projections. Neverthe- less, the projected range reduction in the Caucasus region and northern Iran highlights their importance as long-term refugia, possibly related to higher humidity, stronger environmental and climatic heterogeneity and strong vertical zonation of the forest vegetation. PLOS ONE PLOS ONE | https://doi.org/10.1371/journal.pone.0242280 November 17, 2020 1 / 19 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 OPEN ACCESS Citation: Dagtekin D, Şahan EA, Denk T, Ko¨se N, Dalfes HN (2020) Past, present and future distributions of Oriental beech (Fagus orientalis) under climate change projections. PLoS ONE 15 (11): e0242280. https://doi.org/10.1371/journal. pone.0242280 Editor: Michal Bosela, Technical University in Zvolen, SLOVAKIA Received: June 9, 2020 Accepted: October 29, 2020

Published: November 17, 2020 Copyright: © 2020 Dagtekin et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability Statement: Part of our data is open-source and available to everyone online (http://www.euforgen.org/species/fagus-orientalis/ ), and part of our data belongs to the General Directory of Forestry, Turkey. That is why, we are not allowed to distribute the data ourselves, yet, it is available upon request to the General Directory of Forestry, Turkey. Here is the contact information: Republic of Turkey, General Directorate of Forestry Address: Beştepe Mahallesi So¨ğu¨to¨zu¨ Caddesi No:8/1 06560 Yenimahalle / ANKARA / TURKEY Phone: +90 312 296 40 00

Introduction Climate is changing more rapidly than ever: By 2020, global surface temperature has increased by 1˚C relative to the mean of the years 1951–1980 [1]. Global warming is likely to reach 1.5˚C between 2030 and 2052 if it will be continuing to increase at the current rate [2]. Ongoing cli- mate warming is already causing shifts in species phenology, physiological and behavioural traits, geographical ranges, productivity, and disruption of species interactions [3–5]. Among terrestrial biomes, forests are most affected by ongoing global warming, as it not only impacts the survival rate of tree species, but also forces plants to cope with more frequently occurring extreme events such as severe droughts, floods, wildfires, etc. [6]. Since trees are the dominant species of forest ecosystems, any influence on them would cause cascading effects on the envi- ronment in terms of resource availability, local climate stability and ecosystem services [7]. Any impacts on tree species force dependent organisms to alter their life cycles and even cause their extinction in some cases [8]. Thus, understanding the fate of tree species in response to climate change is of crucial importance for conservation and management practices [4,9]. Broadleaf deciduous forests of Asia Minor occur adjacent to two major biodiversity hotspots, the Mediterranean Basin and the Caucasus [10]. The Euxine-Colchic broadleaf forests of Turkey and Georgia, and the Hyrcanian forests of northern Iran are among the most diverse temperate

forest ecosystems in western Eurasia [11]. High landscape diversity in northern Turkey, Trans- caucasia and northern Iran is buffering the effects of climate change and has created refugia for plant species during past glacial periods [12,13]. However, relatively little is known about how future climates will affect the distribution of broadleaf forests in these regions. Beech (genus Fagus) is widely distributed across Eurasia, where it occurs in environments with sufficient rainfall and mild temperatures [14–17]. One of the most common broadleaf species of western Eurasia is the Oriental beech (Fagus orientalis Lipsky), which plays an important role in forming pure or mixed beech-conifer forests. This species occupies a narrow ecological niche being sensitive to late spring frosts and summer drought [18]. Consequently, it has been identified as especially vulnerable to climate change effects [19–21]. Since F. orien- talis acts as a keystone species in its habitat and serves as an important resource to ecosystem services such as plywood, particleboard, furniture, flooring veneer, mining poles, railway tiles, paper, and firewood [22,23], it is of paramount importance to understand the fate of this spe- cies under climate change. The common approach to investigate the impacts of climate change on plant distributions is using species distribution models (SDMs) [24,25]. These models help to identify regions in an area with changing environmental variables that have similar conditions to localities where the

species has been recorded. By using occurrence (presence/absence) and abiotic environmental data, SDMs are frequently used as a tool for estimating the extent of a species’ range in the future or in the past [25,26]. Although a large number of authors have applied SDM to European woody species extensively in order to understand past and future distributions (e.g. [24–30]), there is a paucity of studies focused on Minor Asian tree taxa (e.g. [31,32]). Likewise, a great number of palaeobotanical studies investigated the late Quaternary history of Fagus sylvatica L. in Europe (e.g. [33–36]), and these studies were integrated with phylogeographic investigations (e.g. [37]). In this study, we investigated the potential distribution of the Oriental beech in the past and in the future. To that end, we applied five SDM algorithms provided with the R-package ‘bio- mod2‘to model potential distributions of this species under the present (1960–1990), past [Last Glacial Maximum (LGM), 21 ka, and mid-Holocene (MH), 6 ka], and future (Represen- tative Concentration Pathway (RCP) 2.6, 4.5 and 8.5 for 2050 and 2070) climatic conditions obtained from MIROC-ESM and CCSM4 global climate models. Specifically, we pursue the following objectives: (i) to model the present distribution of F. orientalis, (ii) to reconstruct PLOS ONE Species distribution modeling of Oriental beech (Fagus orientalis) PLOS ONE | https://doi.org/10.1371/journal.pone.0242280 November 17, 2020 2 / 19 (PBX) Website: https://www.ogm.gov.tr/lang/en/ Pages/Contact.aspx Funding: TD acknowledges funding from the

Swedish Research Council (VR grant no. 2015- 03986). Competing interests: The authors have declared that no competing interests exist.

potential late Quaternary refugia of F. orientalis, (iii) to investigate the late Quaternary history of Fagus based on published palynological data, and finally, (iv) to project the future distribu- tion of F. orientalis under different climate change scenarios. Materials and methods Study region and occurrence data The study region (22˚–54˚ E 35˚–47˚ N) is defined by the observed distribution range of F. orientalis with some extensions for possible future or past range expansions. Occurrence data were obtained from forest management plans of the General Directorate of Forestry of Turkey (GDF) for the distribution of the species in Turkey and the European Forest Genetic Resources Program (F. orientalis–EUFORGEN [38]) for the remaining areas. We merged these two layers of data in ArcGIS (version 10.3.1) and created a single, presence-only dataset (Fig 1) to be used in the models. All occurrence points, obtained from EUFORGEN and GDF, were transformed into 2.5’ spatial resolution using ArcGIS, generating 10,493 presence points for 457,681 raster cells across the study area. Environmental parameters Since the climate is the main abiotic factor shaping species distribution [24], we used nine of 19 bioclimatic variables obtained from WorldClim version 1.4 [39] as the main environmental Fig 1. Modern range of F. orientalis based on occurrence data from EUFORGEN (red) and General Directorate of Forestry of Turkey (green). https://doi.org/10.1371/journal.pone.0242280.g001 PLOS ONE Species distribution modeling of Oriental beech (Fagus orientalis)

PLOS ONE | https://doi.org/10.1371/journal.pone.0242280 November 17, 2020 3 / 19

parameters (see S1 Table). The WorldClim dataset offers global climate model (GCM) simula- tions for the past (LGM and MH) and future (2050 and 2070) climatic conditions. We used MIROC-ESM and CCSM4 global climate simulations since they were the only two GCMs that provided data for all focal periods [40,41]. In addition, for the future simulations, climate sce- narios based on IPCC reports were used [2]. According to these scenarios, in 2100, the expected CO2 levels will be 450 ppm, 650 ppm, and 1350 ppm in optimistic, moderate, and pessimistic scenarios, respectively. In addition, the expected increase in mean annual tempera- ture will be 0.2–1.8˚C, 1.0–2.6˚C, and 2.6–4.8˚C, in optimistic, moderate, and pessimistic sce- narios, respectively [2]. To avoid issues arising from predictor collinearity, we retained only those bioclimatic variables with low correlation (|r| < 0.85) (Table 1, [42,43]). As a result, we used a subset consisting of nine bioclimatic variables (Table 2). In addition, the particular ecol- ogy of beech forests was taken into account when reducing the collinearity within the environ- mental parameters. According to Salamon-Albert et al. [44], temperature variables (annual mean, warmest, and coldest seasonal values) determine the distribution of beech forests. For Oriental beech forests at the latitudes covering our study region both temperature and precipi- tation factors play important roles. In addition, accounting for both temperature and precipi- tation in the model structure is important to consider water and moisture availability [45].

Model structure and analysis SDMs can describe or predict the probability of the presence or absence of a species across environmental gradients or in a specified geographical area [46]. We used the R-package ‘bio- mod2’ version 3.3–7.1 [47] to create potential-distribution maps of F. orientalis. To evaluate model performance, we implemented five different algorithms commonly used for SDMs pro- vided with ‘biomod2’: Generalized linear model (GLM), general additive model (GAM), ran- dom forest (RF), BIOCLIM, and maximum entropy (MaxEnt). MaxEnt uses presence-only data. In contrast, GLM, GAM, and RF require both presence and absence data, and BIOCLIM uses pseudoabsence data [47], thus, we introduced 20,796 random background points to increase precision, as these add uniformity to the model while still subjected to climatic con- straints [45]. Default settings were used for data formatting; 70% of the input data was used as a training sample. The related code is provided in the supplementary (S1 File). Past and future climates were projected through the trained model, which is based on pres- ent conditions with observed distribution data. To evaluate and compare the model perfor- mances we used the Receiver Operating Characteristic Curve (Area Under the Curve, AUC) which takes values between 0 and 1 (AUC > 0.5 indicating that the model performed better than random), higher AUC value meaning better model performance [48]. Model outputs are given as the species’ presence probability in all algorithms.

To display and evaluate possible distributions we used both lowest predicted value [49] and sensitivity–specificity equality [50] methods to set a threshold for the minimum predicted value of observed occurrences, which was calculated as 0.585. We aimed at zero omission by setting this threshold and tried to balance sensitivity and specificity for increased accuracy. Fossil pollen data To compare and test model performances for past projections we compiled published palyno- logical records of Fagus comprising LGM and MH for the range of F. orientalis. Although Fagus pollen had low values in several pollen diagrams pointing to extra-local origin of the mother plant (e.g. [51]), we accepted such records as evidence for presence of oriental beech in the region during target periods (LGM, MH). PLOS ONE Species distribution modeling of Oriental beech (Fagus orientalis) PLOS ONE | https://doi.org/10.1371/journal.pone.0242280 November 17, 2020 4 / 19

Phylogeographic and taxonomic framework and inferred glacial refugia In general, populations of F. orientalis are markedly more differentiated than F. sylvatica L. with FST (coefficient of differentiation among populations) being 0.157 in the former and 0.032 in the latter. In addition, much higher levels of allelic richness in this species indicating that Pleistocene bottlenecks did not deplete the gene pool to the extent they did in the European F. sylvatica [52]. Geographical subgroups of F. orientalis are morphologi- cally and genetically highly distinct [52–56]. Populations ranging from eastern Bulgaria to northern Turkey and the Amanos Mountains form one cluster. Two additional highly dis- tinct clusters involve the Caucasian populations on the one hand and the populations south of the Caspian Sea on the other hand. In contrast, the European F. sylvatica forms a single cluster [52]. In a recent study, Go¨mo¨ry et al. [56] tested speciation scenarios in sub- groups of F. orientalis under an approximate Bayesian framework. The number of genera- tions was used to estimate divergence times between genetically distinct regional populations. Divergence times suggest that European populations of F. sylvatica diverged from Asian Minor populations (F. orientalis) at c. 1.222–0.7 Ma. The Crimean beech is a hybrid between Caucasian F. orientalis and F. sylvatica. It might have originated at around the Eemian interglacial (130–114 ka). Differentiation among the eastern populations hap- pened much earlier: Caspian populations become isolated from Caucasian ones at 2.2–1.6

Ma, and from Turkish ones at 1.9–1.8 Ma. Finally, genetic isolation between populations of F. sylvatica from the Balkans and Central Europe and Apennine occurred much later, at 100 ka and 70 ka. This has important implications for inferring Pleistocene refugia for beech in western Eurasia. Essentially, when considering LGM and MH, in situ refugia must have existed for all geographical subgroups of F. orientalis because these regional groups had been isolated long before these events. Table 1. Correlation matrix between 19 bioclimatic variables. VARIABLES BIO1 BIO10 BIO11 BIO12 BIO13 BIO14 BIO15 BIO16 BIO17 BIO18 BIO19 BIO2 BIO3 BIO4 BIO5 BIO6 BIO7 BIO8 BIO9 BIO1 1,00 BIO10 0,97 1,00 BIO11 0,96 0,87 1,00 BIO12 -0,45 -0,56 -0,29 1,00 BIO13 -0,59 -0,65 -0,46 0,90 1,00 BIO14 -0,21 -0,32 -0,06 0,86 0,60 1,00 BIO15 -0,17 -0,04 -0,31 -0,46 -0,10 -0,77 1,00 BIO16 -0,59 -0,66 -0,46 0,91 0,99 0,60 -0,10 1,00 BIO17 -0,21 -0,32 -0,05 0,87 0,60 0,99 -0,78 0,61 1,00 BIO18 -0,68 -0,68 -0,61 0,82 0,87 0,61 -0,17 0,87 0,61 1,00 BIO19 0,26 0,08 0,44 0,53 0,30 0,64 -0,61 0,30 0,64 -0,02 1,00 BIO2 -0,09 0,00 -0,22 -0,32 -0,25 -0,34 0,50 -0,24 -0,36 -0,23 -0,28 1,00 BIO3 -0,19 -0,29 -0,07 0,23 0,19 0,15 0,05 0,20 0,14 0,01 0,31 0,62 1,00 BIO4 0,19 0,42 -0,09 -0,61 -0,47 -0,55 0,50 -0,47 -0,56 -0,25 -0,64 0,41 -0,44 1,00 BIO5 0,93 0,98 0,80 -0,63 -0,71 -0,40 0,06 -0,71 -0,41 -0,74 0,02 0,20 -0,18 0,49 1,00

BIO6 0,92 0,81 0,99 -0,20 -0,38 0,01 -0,38 -0,38 0,02 -0,54 0,49 -0,37 -0,15 -0,18 0,72 1,00 BIO7 0,03 0,24 -0,23 -0,58 -0,45 -0,56 0,59 -0,45 -0,57 -0,28 -0,62 0,75 -0,05 0,90 0,39 -0,36 1,00 BIO8 -0,15 -0,08 -0,23 0,25 0,33 0,11 0,11 0,32 0,12 0,51 -0,36 -0,08 -0,25 0,25 -0,11 -0,21 0,13 1,00 BIO9 0,81 0,75 0,81 -0,56 -0,67 -0,37 -0,01 -0,66 -0,37 -0,87 0,31 0,14 0,13 0,03 0,78 0,75 0,04 -0,51 1,00 The bioclimatic variables selected for the present study are highlighted in red. https://doi.org/10.1371/journal.pone.0242280.t001 PLOS ONE Species distribution modeling of Oriental beech (Fagus orientalis) PLOS ONE | https://doi.org/10.1371/journal.pone.0242280 November 17, 2020 5 / 19

Results Model evaluation and environmental parameters relevant for Fagus orientalis growth in Asia Minor The results are indicating that models were good in distinguishing presence data from the background based on their AUC values varying from 0.79 to 0.99 (all > 0.7, Table 3, [57]). According to this, the best-fit algorithm was RF and the worst-fit was BIOCLIM. In the follow- ing, we will focus on the RF algorithm provided with ‘biomod2’, which fitted best for F. orien- talis. Successful model results from the present climate could be used to identify and evaluate optimal environmental conditions for F. orientalis growth. The projections for the past, MH and LGM, and for the future, 2050 and 2070, gave different distributions with different climate change scenarios and global climate models (S1–S10 Figs). Past projections Projections for the LGM at 21 ka and the MH at 6 ka from two different global climate models gave consistent results. Past projection outputs of RF are shown with the pollen records col- lected from published palynological studies of Fagus during LGM and MH (Fig 2). The model outputs show that the possible distribution of Fagus was mainly from the mid-Black Sea region to Caucasia and northern Iran to southern Turkmenistan during the LGM (Fig 2A). Addition- ally, RF predicted a refuge area in the Amanos region during the LGM. For the MH, the popu- lations ranging from the mid-Black Sea region to Caucasia remained in these regions, whereas, the Iranian/Turkmen region almost lost its population, except some patches. During the MH,

species distribution shifted towards the west, especially towards Europe and western Anatolia (Fig 2B). Outputs arising from other algorithms are shown in S1–S10 Figs. Late Quaternary palynological records of Fagus in Asia Minor A general feature of the 30 pollen profiles containing Fagus is the sparse occurrence or absence of this genus during the LGM and its more or less continuous presence since about 10 ka (Fig 2 and Table 4). LGM values are usually low indicating long distance dispersal (LDD). One Table 2. Bioclimatic variables (precomputed in WorldClim dataset) used as environmental input in the models. Abbreviations Bioclimatic variables Unit BIO1 Annual Mean Temperature ˚C BIO2 Mean Diurnal Range ˚C BIO3 Isothermality - BIO4 Temperature Seasonality ˚C BIO8 Mean Temperature of Wettest Quarter ˚C BIO9 Mean Temperature of Driest Quarter ˚C BIO12 Annual Precipitation mm BIO15 Precipitation Seasonality mm BIO19 Precipitation of Coldest Quarter mm https://doi.org/10.1371/journal.pone.0242280.t002 Table 3. AUCtest values of all the algorithms, bioclim, maximum entropy (MaxEnt), generalized additive model (GAM), generalized linear model (GLM), and random forest (RF) performed with present climate conditions (1960–1990). Algorithm BIOCLIM MaxEnt GAM GLM RF AUCtest 0.79 0.94 0.96 0.93 0.99 https://doi.org/10.1371/journal.pone.0242280.t003 PLOS ONE Species distribution modeling of Oriental beech (Fagus orientalis) PLOS ONE | https://doi.org/10.1371/journal.pone.0242280 November 17, 2020 6 / 19

exception is the profile from the Sea of Marmara (core MD01-2430) that records beech with moderate abundance during the LGM. No suitable sediments recording the plant cover during the LGM are known from the southern Caspian area. This might be due to the marked drop in sea level (-50 m) of the Caspian Sea during this period [82]. In contrast, MH distribution of F. orientalis was continuous and abundant throughout its present range from eastern Bulgaria, the Sea of Marmara, northern Turkey, Crimea, the western Caucasus, and the southern Cas- pian Sea region. Percentages commonly are way above the proposed threshold values indicat- ing local presence of Fagus reported in Lisitsyna et al. [51]. Moreover, we did not find convincing evidence for the presence of beech during MH outside its present distribution range. Two pollen profiles from Cappadocia and Lake Van have sporadic occurrences of Fagus pollen but these are very few pollen grains and most likely indicate LDD. This is also suggested by the overall composition of these pollen floras (open oak woodland in the Acıgo¨l area of Cappadocia [69], forest steppe in the Lake Van area [70]). Future projections The combined results from both global climate models with different scenarios show that in the future, the main geographical shift of F. orientalis will be towards the northeast of its present distribution. According to the RF model, a severe range contraction of the species is expected in relation to the present distribution provided by the projected model and the

Fig 2. Past range projections with random forest algorithm; here only MIROC-ESM is shown (for CCSM4 see S8 Fig). Occurrence probability of the species is increasing from grey to red (absence to presence) for a) LGM and b) MH. Dots are indicating pollen records (Table 4). https://doi.org/10.1371/journal.pone.0242280.g002 PLOS ONE Species distribution modeling of Oriental beech (Fagus orientalis) PLOS ONE | https://doi.org/10.1371/journal.pone.0242280 November 17, 2020 7 / 19

compiled occurrence data (Table 5 and Fig 1). In addition, a refuge area was projected southwest of the Caspian Sea including the southwestern Caspian coastal area and coastal mountains (Fig 3, for additional algorithms see S1–S10 Figs). The proportion of species presence in the RF MIROC-ESM projections from different scenarios is shown in Table 5. According to the table, the proportion of presence cells is 1.49% in the optimistic, 0.95% in the moderate, and 0.42% in the pessimistic climate scenario. Details of the spatial changes are shown in Fig 3. Table 4. Palaeobotanical records for the LGM and MH [58–81]. Country Locality/site LGM MH Pollen record Fagus presence Bulgaria Black Sea Core GGC-18 No data Present 12 ka Cont. since 11 ka Bulgaria Varna Lake/Core 3 No data Present 8 ka Continuous Bulgaria W Black Sea No data Present 12 ka Increase at c. 7 ka Greece Lesbos/Megali Limni 01 Absent Absent 62–22 ka Present at c. 32 ka, LDD NW Turkey Marmaris/MAR94-5 Present Hiatus 30 ka Continuous NW Turkey Marmaris/Core MD01-2430 Present Present 23 ka Continuous NW Turkey Marmaris/MAR98-12 No data Present c. 17 ka Peak at HM NW Turkey Marmaris/MAR98-13 No data Present c. 18 ka Peak at HM NW Turkey Black Sea/Core B-7 No data Present c. 12 ka Peak at HM NC Turkey Yenic¸ağa, Bolu No data Present c. 12 ka Increase since c. 7 ka NC Turkey Abant Go¨lu¨, Bolu No data Present c. 10 ka Increase since c. 10 ka NC Turkey Black Sea/Core 22-GC3 Present, LDD Present 18 ka Cont., rapid increase at 8.5 ka

NC Turkey Black Sea/Core 22-GC3/8 No data No data 134–119 ka Cont., rapid increase at 126 ka NC Turkey Black Sea/Core 25-GC1 Present No data 63–19 ka No specific information C Turkey Cappadocica/Eski Acıgo¨l I, II No data Present, LDD 15.6 ka Sporadic, LDD SE Turkey So¨ğu¨tlu¨, Van No data Absent 7 ka Sporadic since c. 5 ka, LDD SE Turkey Lake Van/Core 90–04 No data Present, LDD 12 ka Sporadic since 10 ka, LDD Georgia Sukhumi/ Core no 723 No data Present c. 10 ka Continuous Georgia Sukhumi/Dziguta Core 1 Present, LDD No data 48–9 ka Increase at c. 13 ka Georgia NE Black Sea/Core Ak-2575 No data Present 10 ka Continuous Georgia E Black Sea/Gagra, Core no 471 No data Present 10 ka Continuous Georgia Various sites Present No data 18±2 ka No information Iran Lake Urmia Absent Absent 200 ka Absent Iran Caspian Sea/Core CP14 No data Present 5 ka Continuous Iran Caspian Sea/Core GS18 No data Present, LDD 12–2 ka Cont. low, LDD Iran Caspian Sea/Core GS05 No data Present, LDD 15 ka Sporadic since c. 9 ka, LDD NE Iran Kongor Lake No data Present, LDD 6 ka Discont. low, LDD https://doi.org/10.1371/journal.pone.0242280.t004 PLOS ONE Species distribution modeling of Oriental beech (Fagus orientalis) PLOS ONE | https://doi.org/10.1371/journal.pone.0242280 November 17, 2020 8 / 19

Discussion Model evaluation for different SDM algorithms Modelling the spatial distribution of F. orientalis provided details of the species responses to different climate change scenarios and identified possible past refugia by means of projected past distributions and palynological data. A key assumption of SDM with bioclimatic variables is that the ranges of the modelled species are in equilibrium with climate [25, 29, 83]. In gen- eral, the modelled present distributions fitted reasonably well with the observed distribution data, which allowed us to use the trained model for past and future simulations. As mentioned before, all models ran successfully according to their AUC values (Table 3), and the combina- tion of regression and classification approaches, as well as the in-algorithm verification process provided more accurate outcomes in RF and made this algorithm the most successful one [84]. On the other hand, according to the AUC values (Table 3) the least successful algorithm was BIOCLIM. Due to the advanced characteristics and calculations of other algorithms, it was expected that they would give a better fit than BIOCLIM. One problematic algorithm was GAM (S5 and S6 Figs) with a very wide distribution of Oriental beech projected for LGM and MH. Although GAM has the second-best AUC value, we interpreted this as an overfitted dis- tribution. It is known that this algorithm is highly sensitive to large sample size since the fitted functions are not constrained to any functional form when sample size increases [85]. Further,

it is noteworthy that while there are substantial differences between modelling performances from different algorithms, the outcomes are consistent, which is a hint for repeatability [86]. The present distribution of Oriental beech The present distribution of Fagus orientalis in the RF algorithm/MIROC-ESM was largely con- gruent with the observed presence records (Figs 1 and 3A), which is also illustrated by its high AUC value. A range expansion in the European region (Fig 3A) most probably reflects similar present growing conditions of F. orientalis and its sister species F. sylvatica. Hence, it is not surprising that models recognize habitats of F. sylvatica as suitable for F. orientalis [8, 17]. Although there is a close match between observed and modelled distribution, it is interesting to note that the model predicted range expansion to the south of the Pontic Mountains (Kuzey Anadolu Dağları). Under present conditions, this is unlikely as the rain shadow south of the main ridge of the Pontic Mountains creates a sharp change from fully humid Cf climates (see [16, 87]) to Mediterranean Cs and continental, arid BS climates. Whereas, Fagus is replaced by Pinaceae at higher elevations, it is replaced by Quercus species (e.g. Quercus macranthera Fisch. et Mey.) both at higher elevations and in drier areas beyond the main ridge [88]. Table 5. Spatial cell numbers of contraction, overlap, and expansion shown in Fig 3 are quantified in this table. Period Future Scenario Expected Contraction Expected Contraction in model Overlap Expected Expansion Overlap with observed

occurrence data (%) Present – 115 – 10,372 9.444 98.85 2070 RCP 2.6 7,944 7919 3,954 1.558 37.68 RCP 4.5 9,009 8226 2,582 927 24.61 RCP 8.5 9,816 8914 1,087 446 10.36 Overlap of the species occurrence between the observed data and model outputs are given as percentages. In future projections, from optimistic to pessimistic scenario, contractions are increasing, overlap and expansion are decreasing. https://doi.org/10.1371/journal.pone.0242280.t005 PLOS ONE Species distribution modeling of Oriental beech (Fagus orientalis) PLOS ONE | https://doi.org/10.1371/journal.pone.0242280 November 17, 2020 9 / 19

PLOS ONE Species distribution modeling of Oriental beech (Fagus orientalis) PLOS ONE | https://doi.org/10.1371/journal.pone.0242280 November 17, 2020 10 / 19

The situation is similar in the Iranian Alborz Mountains, where Quercus macranthera replaces Fagus at higher elevations [89]. The Amanos Mountains are distinctive, because Ori- ental beech is confined to higher elevations under a distinct microclimate, thus, it is remark- able that the model could identify this relict population here [90]. Past reconstructions and palynological data comparison Probability estimates for occurrences of Oriental beech during the LGM by ‘biomod2’ (Fig 2A) are to some extent congruent with known pollen records. Unambiguous occurrences of beech pollen from LGM deposits are known from the Sea of Marmara, off the central northern Turk- ish Black Sea coast, and from Sukhumi. These occurrences coincide with highly probable model estimates (red in Fig 2A). While estimates on presence in Transcaucasia and southwest of the Caspian Sea are sensible although not supported by palynological data, a refuge area in Turkmenistan at 60˚E east of the Caspian Sea is highly unrealistic in view of the present arid climate in this region. During the LGM, a cold continental climate would not have provided suitable conditions for the growth of beech. ‘biomod2’ reconstructed a further refugium in the Alborz Mountains south of 35˚N. Although no pollen record is available for this region, the Amanos Mountains with a high relief, deep valleys, and slopes facing the sea likely provided a refuge for beech populations during the LGM. The core region of F. orientalis estimated for the MH (Fig 2B) includes the Turkish Black

Sea Coast and the Caucasus. This is also suggested by unambiguous pollen records known from this area (Fig 2B and Table 4). In contrast, ‘biomod2’ with RF predicted a rather limited distribution south of the Caspian Sea, from where abundant pollen records of beech are known as well. Palynological data and paleoclimate modelling [91, 92] suggest that temperate regions during the MH might have been warmer than today during summer and colder during winter because of astronomical forcing. Today, fully humid, warm temperate Cf climates are much less common along the southern Caspian Sea coast than south of the Black Sea and along the northern and southern foothills of the Caucasus ([16]; http://koeppen-geiger.vu- wien.ac.at/kmz/Global_1986–2010_KG_5m.kmz.zip provided by Veterinary University of Vienna, http://koeppen-geiger.vu-wien.ac.at/present.htm). In northern Iran, a stronger Medi- terranean influence is expressed by lower summer precipitation and higher summer tempera- tures (http://koeppen-geiger.vu-wien.ac.at/kmz/Global_1986–2010_KG_5m.kmz.zip). Hence, warmer summers during the Mid Holocene would have increased this Mediterranean influ- ence and resulted in a much-reduced distribution of Oriental beech in the models. As today, these unfavourable conditions might have been compensated by high air moisture and a sky covered with fog for most of the year and particularly during the growing season [89]. How- ever, the predicted presence of beech in the western corner of the Caspian coast makes sense,

as this area–constrained by wind directions and high mountains to the south and a large water body to the north, very likely could have supported growth of beech during the MH (Fig 2B). Future projections Future projections show that F. orientalis will face a dramatic range contraction in the future (Figs 3B–3D and S1–S10). The projections suggest that the species will mainly thrive in the Fig 3. Distribution maps of present and future projections of F. orientalis are shown from the random forests outputs of a) present, b) optimistic future scenario RCP 2.6, c) moderate future scenario RCP 4.5, and d) pessimistic future scenario RCP 8.5. Red areas show the contraction of species occurrence based on observed data and future projections, brown areas show the contraction of species occurrence based on present and future projections, blue areas represent overlap of species occurrence based on observed data and present projection, green areas show the expansion of species occurrence based on observed data and projections. We quantified the cell numbers for each category in Table 5. Spatial resolution: 2.5’. https://doi.org/10.1371/journal.pone.0242280.g003 PLOS ONE Species distribution modeling of Oriental beech (Fagus orientalis) PLOS ONE | https://doi.org/10.1371/journal.pone.0242280 November 17, 2020 11 / 19

Comparison with SDMs of European beech (Fagus sylvatica) Process-based and statistical modelling approaches have been widely used for modelling past and future distribution of a close relative of F. orientalis, Fagus sylvatica. For instance, Gisecke et al. [27] used the process-based bioclimatic model to understand the distribution of F. sylvatica during the MH. In agreement with our findings, they also compared their model results with pollen records and found some mismatches, in particular in the northern and western parts of the distribution. They emphasized the need for testing models against palynological data to increase their predictive value. Likewise, a recent study modelling the modern and past distribution of deciduous oaks (Quercus robur L.) also showed a mismatch between the modelled LGM distribution and the one inferred from pollen data [32]. This study suggested that the more restricted distribution data from the pollen record might be due to competition of deciduous oaks with pines and Mediterranean evergreen oaks. Saltre´ et al. [94] applied a process-based approach coupled with migration models to assess whether the present distribution of F. sylvatica is controlled more by climatic conditions or the migration ability during the MH in Europe. Their findings suggested that the northeast boundary of F. sylvatica distribution is limited by climatic conditions whereas the north- west boundary is limited by migration ability. This is very similar to the situation in F. orien-

talis, where, clearly, the (north)eastern boundary is controlled by climate whereas the west- ern boundary is limited by migration ability. However, in contrast to F. sylvatica, the westward expansion of Oriental beech is limited because it is replaced towards the west by its sister species F. sylvatica. Kramer et al. [95] investigated F. sylvatica distribution under climate change with both process-based and statistical modelling approaches. They found that statistical species distribution models exhibited a better goodness-of-fit for future pro- jections, with the southern limit of the F. sylvatica distribution shifted towards the north. On the other hand, future projections of Dyderski et al. [8] for twelve European forest tree species including European beech suggested that the distribution ranges in western and southern Europe will decrease. This seems plausible in view of the predicted climate change towards summer-dry Cs climates in western and southwestern Europe during 2070–2100 ([96]; Global_Shift_A1FI_1976–2100_30m.kmz). PLOS ONE Species distribution modeling of Oriental beech (Fagus orientalis) PLOS ONE | https://doi.org/10.1371/journal.pone.0242280 November 17, 2020 12 / 19

Conclusions In this study, we investigated the late Quaternary history of F. orientalis utilizing palynological data and past reconstructions and estimated the future distribution of the species in response to global climate change. We compared model outputs to the observed present distribution of F. orientalis and palynological records from LGM and MH. Based on these empirical data, the RF algorithm, among five algorithms in the R-package ‘biomod2’, performed best with more plausible results for both past and future distributions. As in previous studies using F. sylvatica and Quercus robur as focal species, models tend to overestimate occurrence probabilities. For LGM and MH, the results were roughly in accord with palynological data for the Sea of Mar- mara and the Black Sea region, while palynological data did not support potential refugia south of the Caspian Sea during LGM. However, genetically long isolated populations of Ori- ental beech south of the Caspian Sea strongly suggest that they survived the LGM in-situ. For future projections, with continuing climate change, F. orientalis populations in the western Eurasian region will be decreasing with increasing drought; the geographical range will con- tract leaving only two possible refugia, the southern Caspian Sea region (northern Iran) and the Caucasus. This study shows, albeit they can be improved by including other environmental data (e.g. aspect, soil), that our models provide crucial information on the effect of climate change on F.

orientalis in Asia Minor even when only considering the climate. In addition, these results emphasize the importance of combining evidence from modelling with real-world data, such as the rich Quaternary pollen record of Fagus, to understand species-specific responses to changing environments. Hence, these results should be taken into account in order to improve conservation and management plans for this important forest species in western Eurasia. Supporting information S1 Table. 19 bioclimatic variables obtained from WorldClim version 1.4. (PDF) S1 File. R code sheet for biomod2 package to perform SDM. (PDF) S1 Fig. Probability distributions of F. orientalis predicted for the present, MH, LGM, pessi- mistic (RCP 8.5) 2050 and 2070, moderate (RCP 4.5) 2050 and 2070, and optimistic (RCP 2.6) 2050 and 2070 with MaxEnt algorithm and MIROC–ESM GCM. (TIF) S2 Fig. Probability distributions of F. orientalis predicted for the present, MH, LGM, pessi- mistic (RCP 8.5) 2050 and 2070, moderate (RCP 4.5) 2050 and 2070, and optimistic (RCP 2.6) 2050 and 2070 with MaxEnt algorithm and CCSM4 GCM. (TIF) S3 Fig. Probability distributions of F. orientalis predicted for the present, MH, LGM, pessi- mistic (RCP 8.5) 2050 and 2070, moderate (RCP 4.5) 2050 and 2070, and optimistic (RCP 2.6) 2050 and 2070 with GLM algorithm and MIROC–ESM GCM. (TIF) S4 Fig. Probability distributions of F. orientalis predicted for the present, MH, LGM, pessi- mistic (RCP 8.5) 2050 and 2070, moderate (RCP 4.5) 2050 and 2070, and optimistic (RCP 2.6) 2050 and 2070 with GLM algorithm and CCSM4 GCM. (TIF) PLOS ONE

Species distribution modeling of Oriental beech (Fagus orientalis) PLOS ONE | https://doi.org/10.1371/journal.pone.0242280 November 17, 2020 13 / 19

S5 Fig. The probability of F. orientalis predicted for the present, MH, LGM, pessimistic (RCP 8.5) 2050 and 2070, moderate (RCP 4.5) 2050 and 2070, and optimistic (RCP 2.6) 2050 and 2070 with GAM algorithm and MIROC GCM. (TIF) S6 Fig. The probability of F. orientalis predicted for the present, MH, LGM, pessimistic (RCP 8.5) 2050 and 2070, moderate (RCP 4.5) 2050 and 2070, and optimistic (RCP 2.6) 2050 and 2070 with GAM algorithm and CCSM4 GCM. (TIF) S7 Fig. Probability distributions of F. orientalis predicted for the present, MH, LGM, pessi- mistic (RCP 8.5) 2050 and 2070, moderate (RCP 4.5) 2050 and 2070, and optimistic (RCP 2.6) 2050 and 2070 with RF algorithm and MIROC–ESM GCM. (TIF) S8 Fig. Probability distributions of F. orientalis predicted for the present, MH, LGM, pessi- mistic (RCP 8.5) 2050 and 2070, moderate (RCP 4.5) 2050 and 2070, and optimistic (RCP 2.6) 2050 and 2070 with RF algorithm and CCSM4 GCM. (TIF) S9 Fig. Probability distributions of F. orientalis predicted for the present, MH, LGM, pessi- mistic (RCP 8.5) 2050 and 2070, moderate (RCP 4.5) 2050 and 2070, and optimistic (RCP 2.6) 2050 and 2070 with BIOCLIM algorithm and MIROC–ESM GCM. (TIF) S10 Fig. Probability distributions of F. orientalis predicted for the present, MH, LGM, pes- simistic (RCP 8.5) 2050 and 2070, moderate (RCP 4.5) 2050 and 2070, and optimistic (RCP 2.6) 2050 and 2070 with BIOCLIM algorithm and CCSM4 GCM. (TIF) Acknowledgments We are thankful to Abbas Şahin and the General Directorate of Forestry, Turkey, for providing

Fagus distribution data in Turkey. Valuable comments of two anonymous reviewers helped improving the final manuscript. Author Contributions Conceptualization: Dilsad Dagtekin, Evrim A. Şahan, Nesibe Ko¨se, H. Nu¨zhet Dalfes. Data curation: Dilsad Dagtekin, Evrim A. Şahan, Thomas Denk, Nesibe Ko¨se, H. Nu¨zhet Dalfes. Formal analysis: Dilsad Dagtekin, Evrim A. Şahan, H. Nu¨zhet Dalfes. Funding acquisition: Thomas Denk. Investigation: Dilsad Dagtekin, H. Nu¨zhet Dalfes. Methodology: Dilsad Dagtekin, Evrim A. Şahan, Thomas Denk, Nesibe Ko¨se, H. Nu¨zhet Dalfes. Project administration: Nesibe Ko¨se, H. Nu¨zhet Dalfes. Resources: Thomas Denk, Nesibe Ko¨se, H. Nu¨zhet Dalfes. Software: Dilsad Dagtekin, Evrim A. Şahan, H. Nu¨zhet Dalfes. PLOS ONE Species distribution modeling of Oriental beech (Fagus orientalis) PLOS ONE | https://doi.org/10.1371/journal.pone.0242280 November 17, 2020 14 / 19

Supervision: Nesibe Ko¨se, H. Nu¨zhet Dalfes. Validation: Dilsad Dagtekin, Evrim A. Şahan, Thomas Denk, Nesibe Ko¨se, H. Nu¨zhet Dalfes. Visualization: Dilsad Dagtekin, Evrim A. Şahan, H. Nu¨zhet Dalfes. Writing – original draft: Dilsad Dagtekin, Evrim A. Şahan, Thomas Denk, Nesibe Ko¨se, H. Nu¨zhet Dalfes. Writing – review & editing: Dilsad Dagtekin, Evrim A. Şahan, Thomas Denk, Nesibe Ko¨se, H. Nu¨zhet Dalfes. References 1. Vital Signs of the Planet. [cited 07 September 2020]. In: NASA 2020 Global Climate Change [Internet]. Available from: https://climate.nasa.gov/vital-signs/global-temperature/ 2. IPCC, 2018: Summary for Policymakers. In: Global Warming of 1.5˚C. An IPCC Special Report on the impacts of global warming of 1.5˚C above pre-industrial levels and related global greenhouse gas emis- sion pathways, in the context of strengthening the global response to the threat of climate change, sus- tainable development, and efforts to eradicate poverty [Masson-Delmotte V., Zhai P., Po¨ rtner H.-O., Roberts D., Skea J., Shukla P.R., Pirani A., Moufouma-Okia W., Pe´an C., Pidcock R., Connors S., Mat- thews J.B.R., Chen Y., Zhou X., Gomis M.I., Lonnoy E., Maycock T., Tignor M., and Waterfield T. (eds.)]. World Meteorological Organization, Geneva, Switzerland, 32 pp. 3. Walther GR, Post E, Convey P, Menzel A, Parmesan C, Beebee TJ, et al. Ecological responses to recent climate change. Nature. 2002; 416(6879):389–395. https://doi.org/10.1038/416389a PMID: 11919621 4. Parmesan C. Ecological and evolutionary responses to recent climate change. Annu. Rev. Ecol. Evol.

Syst. 2006; 37:637–669. 5. Scheffers BR, De Meester L, Bridge TC, Hoffmann AA, Pandolfi JM, Corlett RT, et al. The broad foot- print of climate change from genes to biomes to people. Science. 2016; 354(6313): 719–730. 6. Lindner M, Maroschek M, Netherer S, Kremer A, Barbati A, Garcia-Gonzalo J, et al. Climate change impacts, adaptive capacity, and vulnerability of European forest ecosystems. For. Ecol. Manage. 2010; 259(4):698–709. 7. Boisvenue C, Running SW. Impacts of climate change on natural forest productivity—Evidence since the middle of the 20th century. Glob. Change Biol. 2006; 12(5):862–82. 8. Dyderski MK, Paź S, Frelich LE, Jagodziński AM. How much does climate change threaten European forest tree species distributions?. Glob. Change Biol. 2018; 24(3):1150–63. https://doi.org/10.1111/gcb. 13925 PMID: 28991410 9. Root TL, Schneider SH. Conservation and climate change: the challenges ahead. Cons. Biol. 2006; 20 (3): 706–8. https://doi.org/10.1111/j.1523-1739.2006.00465.x PMID: 16909561 10. Myers N, Mittermeier RA, Mittermeier CG, Da Fonseca GA, Kent J. Biodiversity hotspots for conserva- tion priorities. Nature. 2000; 403(6772): 853. https://doi.org/10.1038/35002501 PMID: 10706275 11. Temperate broadleaf and mixed forests. [cited March 2019]. In: World Wildlife Fund, WWF [Internet]. Available from: https://www.worldwildlife.org/biomes/temperate-broadleaf-and-mixed-forests 12. Me´dail F, Diadema K. Glacial refugia influence plant diversity patterns in the Mediterranean Basin. J. Biogeogr. 2009; 36(7): 1333–45.

13. Şekercioğlu C¸ H, Anderson S, Akc¸ay E, Bilgin R, Can O¨ E, Semiz G, et al. Turkey’s globally important biodiversity in crisis. Biol Conserv. 2011; 144(12): 2752–69. 14. Yaltırık F, & Elic¸in G. Trees and Shrubs in European Turkey. İstanbul U¨ niversitesi Orman Faku¨ ltesi Der- gisi A. 1982; 32(2). 15. Fang J, Lechowicz MJ. Climatic limits for the present distribution of beech (Fagus L.) species in the world. J. Biogeogr. 2006; 33(10): 1804–19. 16. Kottek M, Grieser J, Beck C, Rudolf B, Rubel F. World map of the Ko¨ppen-Geiger climate classification updated. Meteorol. Z. 2006; 15(3): 259–63. 17. Caudullo G, Welk E, San-Miguel-Ayanz J. Chorological maps for the main European woody species. Data in brief. 2017; 12: 662–666. https://doi.org/10.1016/j.dib.2017.05.007 PMID: 28560272 18. Peters R. Beech Forests: Woody Species Composition, Populations and Spatial Aspects. In: Beech Forests. Dordrecht: Springer;1997. p. 89–130. PLOS ONE Species distribution modeling of Oriental beech (Fagus orientalis) PLOS ONE | https://doi.org/10.1371/journal.pone.0242280 November 17, 2020 15 / 19

19. Ko¨ se N, Gu¨ner HT. The effect of temperature and precipitation on the intra-annual radial growth of Fagus orientalis Lipsky in Artvin, Turkey. Turk. J. Agric. For. 2012; 36(4): 501–9. 20. Haghshenas M, Mohadjer MR, Attarod P, Pourtahmasi K, Feldhaus J, Sadeghi SM. Climate effect on tree-ring widths of Fagus orientalis in the Caspian forests, northern Iran. Forest Sci. Technol. 2016; 12 (4): 176–82. 21. Martin-Benito D, Pederson N, Ko¨ se N, Doğan M, Bugmann H, Mosulishvili M, et al. Pervasive effects of drought on tree growth across a wide climatic gradient in the temperate forests of the Caucasus. Glob. Ecol. Biogeogr. 2018; 27(11): 1314–25. 22. Kandemir G, Kaya Z. EUFORGEN [European Forest Genetic Resources Programme] technical guide- lines for genetic conservation and use for oriental beech (Fagus orientalis). EUFORGEN Technical Guidelines for Genetic Conservation and Use. 2009. 23. Bozkurt AY, Erdin N. Ağac¸ Teknolojisi Ders Kitabı, İstanbul U¨ niversitesi Orman Fak. Orman End. Mu¨ h. Bo¨ l., 372s. İstanbul. 1997. 24. Pearson RG, Dawson TP. Predicting the impacts of climate change on the distribution of species: are bioclimate envelope models useful? Glob. Ecol. Biogeogr. 2003; 12(5): 361–71. 25. Pearson RG. Species’ distribution modeling for conservation educators and practitioners. Synthesis. Am. Mus. Nat. Hist. 2007; 50: 54–89. 26. Elith J, Leathwick JR. Species distribution models: ecological explanation and prediction across space and time. Annu. Rev. Ecol. Evol. Syst. 2009; 40: 677–97.