Journal of Plant Ecology VolumE 9, NumbEr 4, PagEs 421–433 august 2016 doi:10.1093/jpe/rtv069 advance access publication 14 october 2015 available online at www.jpe.oxfordjournals.org transition along gradient from warm to mesic temperate forests evaluated by gamm Andraž Čarni1,2,3,*, Vlado Matevski3,4, Nina Juvan5, Mitko Kostadinovski4, Petra Košir6, Aleksander Marinšek7, Andrej Paušič8 and Urban Šilc1,9 1 Abstract Aims the aim of the study was to discover what set of variables best explains the transition from warm to mesic forest vegetation. based on various variables grouped into sets (geomorphological, ecological, structural, soil characteristics and chorological), six models were built and tested by generalized additive mixed models (gamms). We assumed that each set of variables has different explanatory power. our aim was to compare the six different models (sets of variables), to test which model best explains the species turnover in forest communities along the transition between warm and mesic temperate forests and to try to find reasons for the different explanatory power of the models. Methods the research took place in the southern part of the balkan Peninsula. Field sampling was done according to standard methods. the gradient from warm to mesic forests was defined as the turnover of species and evaluated by projection of samples on the first unconstrained DCa axis. geomorphological, ecological, structural and INtroDuCtIoN The species composition of a biological community results from the interplay between ecological and evolutionary soil characteristics, together with chorological sets of variables, were regressed on the turnover of species composition. based on the five sets of variables, six models were constructed and tested by generalized additive mixed models. Important Findings Ecological conditions best explain the change of forest communities along the gradient; evolution and the development of vegetation reflected in chorotypes are also of high importance; geomorphology and structure seem not to change so dramatically and soil shows the least significant differences of all. Ecological variables are the most important set of variables in the transition between warm and mesic temperate forests but eco-evolutionary dynamics after the Pleistocene should also be taken into consideration. Keywords: balkan, chorotypes, life forms, ecology plant ecology, soil science, vegetation, evolution, refugium received: 16 april 2015, revised: 18 september 2015, accepted: 8 october 2015 processes (Ricklefs 1987).The ecological conditions of a site are reflected in various sets of variables, such as geomorphology, soil measurements, bioindicator values and plant functional traits. The last named is defined as non-phylogenetic © The Author 2015. Published by Oxford University Press on behalf of the Institute of Botany, Chinese Academy of Sciences and the Botanical Society of China. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Downloaded from https://academic.oup.com/jpe/article/9/4/421/2222523 by guest on 03 June 2022 Institute of Biology, Research Centre of the Slovenian Academy of Sciences and Arts, Novi trg 2, SI-1000 Ljubljana, Slovenia University of Nova Gorica, Vipavska 13, SI-5000 Nova Gorica, Slovenia 3 Macedonian Academy of Sciences and Arts, Bul. Krste Misirkov, 2, P.O. Box 428, MK-1000 Skopje, Republic of Macedonia 4 Faculty of Natural Sciences and Mathematics, Institute of Biology, University of Ss. Cyril and Methodius, Gazi Baba b/b - P.O. Box 162, MK-1000 Skopje, Republic of Macedonia 5 Anton Melik Geographical Institute, Research Centre of the Slovenian Academy of Sciences and Arts, Novi trg 2, SI-1000 Ljubljana, Slovenia 6 Faculty of Mathematics, Natural Sciences and Information Technologies, University of Primorska, Glagoljaška 8, SI-6000 Koper, Slovenia 7 Slovenian Forestry Institute, Večna pot 2, SI-Ljubljana, Slovenia 8 Kranjčeva 8, SI-9220 Lendava, Slovenia 9 Biotechnical Centre Naklo, Strahinj 99, SI-4202 Naklo, Slovenia *Correspondence address. Institute of Biology, Research Centre of the Slovenian Academy of Sciences and Arts, Novi trg 2, SI-1000 Ljubljana, Slovenia. Tel: +386-1-470-63-11; Fax: +386-1-425-77-97; E-mail: carni@zrc-sazu.si 2 422 are non-linear (Danz et al. 2013). Our analysis was performed by generalized additive mixed models (GAMMs), which allow combined analysis of both linear and non-linear variables in a single model. Modelling enables a set of variables to be combined in a single model and then various models to be compared with each other (Paušič and Čarni 2013; Seddon et al. 2015). GAMMs (Wood 2006) were then used to explain the turnover of species along projection of samples on the first unconstrained DCA axis as a proxy for turnover of species composition along a gradient from warm to mesic temperate forests (AX1). Two factors contributed to the selection of this approach. First, inspection of the raw data showed that non-linear models were appropriate. Generalized additive models (GAMs) are data-driven rather than model-driven. These non-parametric regression models allow determination of the response curve shape from the data, instead of fitting an a priori parametric model that is limited in its available shape of response (Armas et al. 2011; Caprio et al. 2009; Wood and Augustin 2002). GAMs use smoothing functions to fit non-linear response curves (Haslem et al. 2011; Wood 2006). Secondly, the clustered distribution of relevés suggested a mixed model approach. Mixed models are recommended when data are structured by some factor that introduces systematic variation of the potential influence over the relationship between predictor and response variables (Haslem et al. 2011; Zuur et al. 2009). Models were applied using the R programme and mgcv package (R Development Core Team 2013, Vienna, AT; Wood 2006). Based on five sets of variables (geomorphological, ecological, structural, soil characteristics and chorological) and model containing all variables, six models were built and tested by GAMMs. We assumed that each set of variables has different explanatory power. Our aim was to compare different models (sets of variables), to test which model best explains the species turnover in forest communities along the transition between warm and mesic temperate forests and to try to find reasons for the different explanatory powers of the models. matErIals aND mEtHoDs Study area Figure 1: geographical location of the study area in the central part of the Balkans: the Galičica mountain range. We chose a study area of about 250 km2 in the central part of the Balkan Peninsula: the Galičica mountain range, situated between Lake Ohrid (altitude of 695 m) and Lake Prespa (850 m), with a highest point of 2255 m (Fig. 1). The mountain ridge extends in a south–north direction, so the influence of aspect is reduced. Slopes to the lakes are steep, while the upper part is built as a karstic plateau. Since the bedrock is limestone, higher biodiversity resulting from evolutionary and historical processes can be expected (Chytrý et al. 2003; Ewald 2003a). The climate is moderate continental, with 780 mm of precipitation and an average temperature of 11.2°C. The climate above an altitude of 1100 m can be treated as cool continental, with 1067 mm of precipitation Downloaded from https://academic.oup.com/jpe/article/9/4/421/2222523 by guest on 03 June 2022 groups of plant species sharing a similar function at an organismic level and having the same response to environmental factors (Cornelissen et al. 2003; Kraft et al. 2015; Pignatti et al. 2005). At the same time, the importance of the species pool, which reflects the evolutionary and historic angle, must also be taken into consideration (Cubizolle et al. 2014). There are two main phylogeographic patterns (Feliner 2014) that influence the appearance of mesic and warm forests of the region. The first is climatically driven north–south shifts by which species responded to climatic oscillation during the Pleistocene and afterwards. If the area is mountainous, latitudinal and altitudinal interactions should be included (Kropf et al. 2008). The other determinant of the present vegetation is the appearance of potential multiple refugia in the region (Frajman et al. 2007; Kučera et al. 2010; Surina et al. 2011). The research area lies in the transition zone between the Mediterranean and Circumboreal floristic regions (Quézel and Médail 2003) (Fig. 1). The influence of the Mediterranean climate, i.e. mild rainy winters and hot dry summers, and of the Continental one with cold, snowy winters and hot and dry summers, can be recognized in the area. Contact between two large groups of forests can be found, which forms zonal vegetation in the area, as warm (dominated by Quercus sp. div., Ostrya carpinifolia, Carpinus orientalis, Juniperus excelsa) and mesic (dominated by Fagus sylvatica, Carpinus betulus, Acer obtusatum) deciduous forests. Study of this contact is especially important in the light of the foreseen climatic changes, which will induce further changes in the vegetation cover of the region (Čarni and Matevski 2015; Ruiz-Labourdette et al. 2012). We tried to evaluate relationships among geomorphological, ecological, structural, chorological and soil characteristics. It would be possible to study the data by linear regression, which is one of the most frequently used techniques for developing prognostic models (Lang and Ewald 2014; Reger et al. 2011; Wehn et al. 2014). Although this method is exceptionally accurate, it is practically impossible to develop models by means of linear regression when the observed variables Journal of Plant Ecology Čarni et al. | Warm and mesic forests and an average temperature of 7.4°C. Summer droughts appear in the lowlands but they are mitigated at higher altitudes, although winter frost is more severe at higher altitudes (Acevski 2001; Filipovski et al. 1994; Filipovski 1995; Lazarevski 1993). Vegetation sampling (relevés) We measured the geomorphological variables in the field: altitude in meters above sea level, aspect is direction of the slope and is nominal variable, slope is a slope of terrain (in percentage), bare rock is percentage bedrock on the surface not covered by soil or vegetation. Since aspect is not a quantitative variable, we transformed aspect into angle to a northern direction and the cosine function of this angle was used in the analysis. Ecological data The second model (Model 2) was composed of ecological variables. In order to evaluate other ecological factors (temperature, moisture, light, continentality, soil reaction and nutrients), Pignatti indicator values (Pignatti et al. 2005) were calculated using the JUICE program (Tichý 2002). We calculated them on basis of presence–absence data (unweighted). They were determined for Italian flora but because historical and ecological similarity they can also be used for species on Balkan Peninsula (Bergmeier and Dimopoulos 2008; Košir et al. 2008; Marinšek et al. 2013). Although some values for some local species (e.g. Festuca galicicae, Abies borisii-regis) are missing, this does not greatly influence the results (Ewald 2003b). Turnover of species in forest communities along the gradient The DCA (Fig. 2) was performed in the CANOCO version 4.5 program with downweighting of rare species (ter Braak and Šmilauer 2002). The eigenvalues of the first two axis are 0.693 and 0.476. We can detect along the first axis the gradient from mesic temperate forest to warm (thermophilous) forests. On the other axis, the mesic temperate forests show a rather unique groups, but we can detect the separation between the warm (thermophilous) forest thriving on steep slopes on relatively higher altitudes (Quercus pubescens and Ostrya carpinifolia dominated) and lowland (Carpinus orientalis and Juniperus excela dominated) warm (thermophilous) forests. The multivariate data of species composition is simplified by using first (unconstrained) ordination axis as continuous proxy of the gradient from mesic temperate forest to warm (thermophilous) temperate forest; the scores of samples on the first axis are therefore used as a proxy for species turnover along the gradient. Geomorphological data The first model (Model 1) was composed of four geomorphological variables (altitude, aspect, slope and bare rock). Figure 2: ordination diagram showing the result of DCA analysis of communities. Eigenvalues of the first two axes are 0.693 and 0.476, respectively. Downloaded from https://academic.oup.com/jpe/article/9/4/421/2222523 by guest on 03 June 2022 Prior to field sampling, unsupervised classification of satellite imaging of the area was performed (Aronof 2005; Halder et al. 2011). We thus obtained basic information about the diversity of forest vegetation, because we planned to sample all forests types. The field sampling was done according to the standard method for vegetation elaboration (Braun-Blanquet 1964) during two sequential years (2009 and 2010); sampling plots with an area of 100 m2 were located in the centres of larger polygons of forest communities obtained by satellite image classification. In some cases, the satellite image classification produced single polygons dominated by several forest types or the forests were cut during the recent period. Those plots were neglected and we chose (homogenous) plots dominated by single well developed forest type. We took up to 10 samples, most often five, of each forest type determined by classification of satellite images. On each sampling plot, we made a floristic inventory of all presented plant species and estimated their abundance and cover according to seven grade scale according to Braun-Blanquet scale. Such samples are termed relevés (Braun-Blanquet 1964). A total of 105 relevés were made. Relevés were stored in the TURBOVEG database (Hennekens and Schaminée 2001). Raw data are published in a book presenting the forest vegetation of the region (Matevski et al. 2011). 423 424 Structural and functional data This model (Model 3) reflects the visual (number of species, structural and functional) appearance of forests. Number of species is the sum of species in each sampling plot (relevé). Structure is represented by percentage of cover of tree, shrub, herb and moss layers and was visually estimated on each sample plot. Functional data are represented by plant traits (phanerophytes, nano-phanerophytes, chamaephytes, hemicryptophytes, geophytes and therophytes) and derived from species composition and reflect ecosystem functioning (Campetella et al. 2011; Raunkiaer 1934). Soil data Chorotypes We also determined the chorological spectrum of relevés and integrated them into the fifth model (Model 5). Since many classification schemes exist, we accepted the proposed scheme of Pignatti (1982). This scheme is based on actual distribution of plants and was prepared for Italian flora. As this scheme does not encompass the species distributed only in the Balkan Peninsula, we added those species (Gajić 1980) and named them Balkan species. Species were thus divided into the following categories: endemic, steno-mediterranean, euri-mediterranean, mediterranean-montane, eurasiatic, atlantic, orophytes of south Europe, boreal, widely distributed and Balkan species. The species were divided into chorotypes according to various floras (Hegi 1946‒2003; Jordanov 1963‒1982; Josifović 1970‒1977; Diklić 1986; Micevski 1985‒2005; Strid and Tan 1997, 2002). All variables We have taken all variables used in all partial models and built the last model (Model 6). Statistical analysis and modelling techniques Bivariate correlations between different predictor variables were calculated by applying the Spearman rank-correlation test in STATISTICA (StatSoft, inc. 2011). In order to avoid colinearity among the explanatory variables in the model, only one of the highly correlated variables was selected (Alahuhta et al. 2011; Paušič and Čarni 2013). The correlation coefficient was set at ρ > |0.5| to consider the variables, except for the more correlated ecological variables for which the coefficient was set at ρ > |0.8|. Such a correlation coefficient was chosen to include as many variables as possible in a model. All sets of variables included a total of 30 of the potential 44 variables (Table1). In the geomorphological set, all four variables were included: altitude, aspect, bare rock and slope; in the ecological set, four variables, temperature, continentality, soil reaction and nutrient, were included and two, light and moisture, were excluded; in the structural set, six variables were included, chamaephytes, cover of moss layer, cover of tree layer, nano-phanerophytes, number of species, phanerophytes and therophytes and four, cover of shrub layer, cover of herb layer, geophytes and hemicriptophytes, were excluded; in the soil set, seven variables, pH, organic matter, C/N, clay, Mg, K and Na, were included and six, C, N, Ca, H, silt and sand, were excluded; in the chorotype set, eight variables, atlantic, Balkan, boreal, Eurasian, mediterranean-montane, south European oreophytes, steno-mediterranean and wide amplitude, were included and two, endemic and eurimediterranean, were excluded. Descriptive statistics were calculated for each variable included in the models (minimum, maximum, mean and standard deviation) (Table 1). The variables were selected using stepwise regression, which is a combination of backward elimination, whereby explanatory variables are removed from the model, and forward selection, whereby explanatory variables are added to the model. The stepwise regression stops when maximized residual deviance is achieved and new variables do not contribute to the residual deviance (Alahuhta et al. 2011; Faraway 2005). Analysis of each individual variable included in a specific model was carried out. We observed the response of variables in a model (linear–non-linear) and the logic of each variable’s response with regard to the relative measured frequency of the independent variable (AX1). While there is no universal measure of model performance, the use of multiple measurements is common and more objective (Aertsen et al. 2010; Dawson et al. 2007). One of the most powerful approaches for model selection (Burnham and Anderson 2002) is Akaike’s information criterion (AIC; Akaike 1973). The model with the lowest AIC best explains the data (Jiménez-Alfaro et al. 2014b). The Bayesian information criterion (BIC) is a similar measure, which was also tested. Another criterion of model performance that was used is the adjusted coefficient of determination (adjusted R2; adj R2). It is a modification of R2 (Pearson 1896), which also adjusts for the number of explanatory terms used in a model. Unlike R2, adjusted R2 increases only if the new term improves the model more than would be expected by chance. Adjusted R2 can be negative and will always be less than or equal to R2 (Aertsen et al. 2010). All six models were validated using the percentage of explained deviance (D2; Alahuhta et al. 2011). For the occurrence data, evaluation was also based on the area under Downloaded from https://academic.oup.com/jpe/article/9/4/421/2222523 by guest on 03 June 2022 Since the ecological set of variables is derived from species composition, we distinguished them from soil samples that were measured in the field. The fourth model (Model 4) was composed of soil samples collected at each sample plot at a depth from 0 to 10 cm to give additional information about the soil characteristics. The samples were treated in accordance with ISO 11464 and the following parameters were measured pH (CaCl2) (ISO 10390), organic carbon was determined in accordance with the Walkley–Black method (ISO 14235); total nitrogen (N) content according to modified Kjeldhal method (ISO 11261); exchangeable basic cations (Ca2+, Mg2+, K+, Na+) were analyzed using atomic absorption spectrophotometry. Journal of Plant Ecology Linear R2 F P value 0.627 119.5 <2e−16*** 0.6198 Model 1 Non-linear 0.28 18.23 4.9e−8*** 0.28 Model 1 Decrease 0.0156 2676 0.105 0.01364 21.73 Increase 0.16 11.09 0.00162** 0.1596 Model 1 6.14 0.72 Increase 0.0156 2,679 0.105 0.9003 Model 2 4.97 0.29 Non-linear 0.69 101.6 <2e−16*** 0.668 Model 2 5.06 4.2 0.57 Non-linear 0.103 7683 0.00116** 0.0724 Model 2 4.42 6.44 4.42 0.94 Non-linear 0.758 124.9 <2e−16*** 0.7398 Model 2 Value 10 64 36.49 12.37 Non-linear 0.281 41.98 3.22e−9** 0.2806 Model 3 % 0 95 63.43 25.96 Non-linear 0.299 45.74 8.36e−10*** 0.2987 Model 3 Moss_layer % 1 80 14.86 15.47 Non-linear 0.223 59.9 0.000722*** 0.2229 Chamaephytes Cham % 1.64 21.95 8.14 4.58 Increase Nano-Phanerophytes N_Phan % 1.56 20 6.11 4.11 Non-linear Phanerophytes Phan % 8.7 69.57 33.41 11.2 Therophytes Thero % 1.79 46.15 9.58 9.84 pH pH Value 4.6 8 6.17 Organic matter Org_mat % 1.5 62.7 18.81 C/N CN % 1.3 25.9 Clay Clay % 4.3 57.3 Mg Mg_per % 0.3 K K_per % Na Na_per Stenomediterranean Mediterraneanmontane Variable Abbreviation Unit Minimum Maximum Mean SD Trend Altitude Altit Meter 719 1810 1197.65 280.43 Non-linear Aspect Asp Cos −1 1 0.16 0.67 Slope Slope Degrees 2 43 18.95 11.71 Bare rock Bare_rock % 1 80 24.43 Temperature Temp Value 4.65 7.59 Continentality Cont Value 4.3 5.67 Soil reaction SoilR Value 2.82 Nutrient Nutr value Number of species NoSpe Tree layer Tree_layer Cover moss layer GAMM R2 Model selection Geomorphology Čarni et al. Table 1: descriptive statistics of the variables included in the model selection | Warm and mesic forests Ecology Structure 0.0385 4657 0.0336* −0.0105 0.192 0.663 0.03851 Model 3 Increase 0.0562 7268 0.0082** 0.05624 Model 3 Increase 0.194 18.35 5076e−5*** 0.194 Model 3 0.74 Increase 0.268 11.61 6.53e−7*** 0.1903 Model 4 12.22 Non-linear 0.043 1,916 0.151 −0.00979 14.1 3.6 Decrease −0.00978 0.002 0.962 −0.00978 29.84 10.25 Increase −0.00895 0.123 0.726 −0.00895 Model 4 34.1 5.51 4.39 Increase 0.0712 8,988 0.00341** 0.07124 Model 4 0.6 3.5 1.7 0.66 Non-linear 0.32 15.77 1.64e−8*** 0.2199 Model 4 % 17 72 48.02 10.07 Non-linear 0.758 124.9 <2e−16*** 0.010162 Model 4 SM % 1.79 22.73 6.38 3.91 Non-linear 0.116 12.12 0.000799*** ME.MO % 1.56 17.07 5.45 3.23 Decrease 0.0348 0.695 0.475 Eurasian EURIAS % 23.53 88.24 53.23 13.07 Decrease 0.298 45.56 Atlantic ATLANT % 1.56 6.67 2.83 1.07 Non-linear 0.0564 2599 South European orophytes ORO.S.EU % 1.75 20 6.01 4.28 Non-linear 0.477 Boreal BORE % 1.64 26.67 9.36 4.65 Decrease Wide amplitude AMP % 0 12.5 2.51 2.81 Non-linear Balkan BALK % 1.85 30.77 5.91 5.7 Non-linear −0.0105 Soil Chorotypes 0.1155 Model 5 −0.01058 Model 5 8.92e−10*** 0.2979 Model 5 0.12 0.05641 Model 5 61.29 6.46e−11*** 0.4773 Model 5 0.133 16.26 0.00011*** 0.1355 Model 5 0.0835 6578 0.0129* 0.08352 Model 5 0.262 7569 0.00607** 0.2357 Model 5 425 The trends of attributes were deduced from the scatterplots (see supplementary Appendix 1, for examples Fig. 4 in the text). Significance codes: ‘***’P < 0.001; ‘**’P < 0.01; ‘*’P < 0.05. Downloaded from https://academic.oup.com/jpe/article/9/4/421/2222523 by guest on 03 June 2022 426 rEsults Considering floristic composition, the DCA diagram (Fig. 2) shows the arrangement of the relevés in two-dimensional space, along two most significant axis. Along the first axis (AX1), a clear zonation of forests from the submediterranean influenced warm temperate group dominated by Juniperus excelsa, on the one hand, through more continental warm temperate forests dominated by various Quercus species and oriental hornbeam (Carpinus orientalis), to mesic temperate forests dominated by beech (Fagus sylvatica), hornbeam (Carpinus betulus) and ravine forests on the other. This allows the scores of samples on the first axis to be considered as a proxy for a transitional gradient from warm to mesic temperate forests. As variables present several complex phenomena, we grouped them into five sets of variables (i.e. geomorphology, ecology, structure, soil, chorotypes) and, after excluding highly correlated variables in each set, six models were built. The details are presented in Tables 1 and 2, Fig. 3 and supplementary Appendix 1. The first model (Model 1 – Geomorphology) comprises geomorphological properties. The second model (Model 2 – Ecology) comprises average ecological indicator values. The third model (Model 3 – Structure) encompasses the structural and functional characteristics. The fourth model (Model 4 – Soil) comprises soil properties, while the fifth model (Model 5 – Chorotypes) encompasses the origin of species found on sample plots. At last, we joined all variables into one model (Model 6 – All variables) (Table 1). The various explanatory powers for each model were obtained through the inclusion and exclusion of individual variables from each model (supplementary Appendix 2). In each model, we therefore decided on the combination of variables that gives the best results. The selected variables and the models’ explanatory power are presented in Table 2. As seen in Table 2, models differ considerably with regard to their explanatory power. Model 6 (Pearson’s correlation coefficient is 0.99), in which all variables were used, has the highest predictive power. Predicted AX1 values in this model best fit the actual AX1 data (Fig. 4). This model also has the highest percentage of explained deviance (D2) (Table 2). The next highest explanatory power has Model 2 (0.98) in which ecologycal data were used. Model 5 (0.94), in which chorotypes were used, and Model 1 (0.89), with geomorphology data, come close to Model 2 (Table 2). Model 3 (0.85), which comprises structure and function, and Model 4 (0.77) based on soil parameters have the poorest explanatory power and are therefore the least accurate explanatory model for the evaluation of AX1. Estimation of the models’ adequacy (their explanatory power) is fairly abstract information, calculated with selected measured statistical coefficients that (each in turn) describe the power of an individual model. Therefore, there is often a dilemma as to whether the result is truly satisfactory. In order to provide a more illustrative presentation of the models’ explanatory power, a scatter diagram was created. The original input data were entered into each of the six models, presented on the x axis of the diagram and the result of prediction of each model is presented on the y axis. For the sake of clarity, only the trend of the input function was preserved for each model, while the values for relevés were excluded (Fig. 4). DIsCussIoN The ordination diagram (Fig. 2) shows on the first axis a clear distinction between warm and mesic temperate forest. Many tree species can be found as dominants within warm temperate forests, such as various Quercus species: Quercus cerris, Q. frainetto, Q. petraea, Q. pubescens and Q. trojana, as well as Ostrya carpinifolia and Carpinus orientalis. Juniperus excelsa can also be added to this group. This area is the richest region in tree species in Europe (Adams 2009). The second DCA axis reflects higher diversity of warm temperate forest then mesic ones. Besides the climate the reason is that there is one of the main Pleistocene refugial areas of warm (thermophilous) forest vegetation in Europe (Médail and Diadema 2009; Petit et al. 2002, 2008). Here, we can find origin of many thermophilous species and species that could potentially shift to the north under foreseen climate change (Čarni and Matevski 2015). The mesic temperate forests are mainly dominated by Fagus sylvatica, occasionally co-dominated by Abies borisii-regis. The main refugial area of beech forest during Downloaded from https://academic.oup.com/jpe/article/9/4/421/2222523 by guest on 03 June 2022 the curve (AUC) in receiver operating characteristic (ROC) plots. Model accuracy based on AUC was calculated with the R package pROC. (Alahuhta et al. 2011; Swets 1988). Additionally, visual validation of the response curves and spatial predictions were conducted according to standard procedure (Lehmann et al. 2002; Saatkamp et al. 2010; Tamstorf et al. 2007) (Fig. 4). The explanatory power of the models was tested with the input of the dependent variables that explain AX1. The command predict in the vegan package (Oksanen et al. 2013) in R package was used. Real data for variables were entered into the models. The result was the AX1 values calculated by each of the six models. Pearson’s correlation coefficient (r) was used to compare the accuracy of the models. The correlations between the real and the calculated AX1 were assessed in STATISTICA. Pearson’s correlation coefficient is represented in the same way as a correlation coefficient used in linear regression, ranging from −1 to +1. A value of +1 is the result of a perfect positive relationship between two or more variables (Nie et al. 2011). The explanatory power of an individual model was presented in a scatter diagram that incorporated real data of AX1 and the results calculated by the models (Fig. 4). Journal of Plant Ecology Čarni et al. | Warm and mesic forests 427 Table 2: generalized additive mixed models (GAMMs) describing differences in predictions for AX1 through different explanatory variables Model and terms Estimate F value t value P value Model 1: Geomorphology Bare_rock 0.020234 5.988 199.319 <2e−16*** s(Asp) 7.163 0.00215** Model 2: Ecology 1.11857 14.448 9.841 s(SoilR) 5.349 s(Nutr) 66.017 r 221.6461 240.1569 0.766 100 78.8 0.89 0.952 100 95.9 0.98 0.022759* 0.0103* Phan 0.012321 1.525 0.1304 Thero 0.059414 6.226 1.21e−08*** s(NoSpe) 24.85 6.89e−11*** s(Tree_layer) 45.97 5.48e−10*** Model 4: Soil 0.587198 3.960 −0.003517 −0.374 0.060097 2.642 20.135 s(Na_per) 5.119 100 73 0.85 290.42 313.8665 0.546 100 58.8 0.77 209.052 248.8614 0.84 100 87.9 0.94 58.748 143.3686 0.967 100 98.2 0.99 0.709157 0.009704** 4.94e−09*** −0.043831 −2.507 EURIAS −0.084074 −16.641 BORE −0.087124 −7.106 8.203 0.0139* <2e−16*** 2.57e−10*** 0.00516** 4.314 0.01015* s(ORO.S.EU) 29.375 1.14e−12*** s(AMP) 10.252 0.00186** s(BALK) 27.247 9.92e−12*** Model 6: All variables 0.0053958 0.685 0.0059** ME.MO s(ATLANT) 283.4433 0.000148*** Model 5: Chorotypes s(SM) 259.5577 2.276 0.025440* 1.99e−09*** Temp 0.8157872 6.738 Cham −0.0020391 −0.250 0.803125 Phan −0.0020771 −0.701 0.485578 Thero 0.0140625 3.468 0.000834*** CN 0.0001518 0.019 0.984933 ME.MO 0.0013124 0.126 0.899678 EURIAS −0.0040922 −1.038 0.302346 s(Altit) 12.380 0.000703*** s(Cont) 14.211 0.000299*** s(SoilR) 0.268 s(Nutr) 13.408 0.605857 0.000433*** s(NoSpe) 4.506 0.036700* s(Tree_layer) 0.546 0.461877 s(Moss_layer) 5.948 0.002965** Downloaded from https://academic.oup.com/jpe/article/9/4/421/2222523 by guest on 03 June 2022 2.616 s(K_per) 85.86176 2.94e−13*** 0.042584 Bare_rock D2 (%) <2e−16*** Cham Mg_per AUC (%) 0.000107*** Model 3: Structure Clay Adj R2 61.97611 s(Cont) pH BIC 3.4e−08*** s(Altit) Temp AIC 428 Journal of Plant Ecology Table 2: Continued Model and terms Estimate F value t value P value s(N_Phan) 8.416 0.004737** s(Org_mat) 2.090 0.152015 s(SM) 6.700 0.011347* s(ORO.S.EU) 5.061 AIC BIC Adj R2 AUC (%) D2 (%) r Non-bracketed predictors were fitted as typical parametric terms; s(…) denotes predictors fitted as non-parametric smoothing terms. GAMMs use F values to test the significance of non-parametric smoothing terms, while t values are used for parametric terms. Differences in explanatory power based on calculated Akaike’s information criterion (AIC), Bayesian information criterion (BIC), adjusted coefficient of determination (adj R2), area under the curve (AUC) of a receiver operating characteristic plot (ROC), percentage of explained deviance (D2) and Pearson’s correlation coefficient (r). ‘***’P < 0.001; ‘**’P < 0.01; ‘*’P < 0.05. Downloaded from https://academic.oup.com/jpe/article/9/4/421/2222523 by guest on 03 June 2022 Figure 3: example of linear or non-linear response curves between explanatory variable and AX1. The most significant predictors in each model are presented. Around the smooth terms, their 95% confidence intervals are shown. Each short bar on the x axis indicates the distribution of a sample site along that variable. Pleistocene was in the north-western Balkan and the refugia in the southern Europe did not contribute to the population of the continent (Brus 2010; Magri et al. 2006). We consider that Carpinus betulus dominated forests found in humid valleys and Acer obtusatum and Corylus colurna dominated forests in ravines also belong to this group (Košir et al. 2008, 2013). The first model is based on a geomorphological set of variables (Tables 1 and 2, Fig. 3, supplementary Appendix 1). It can be seen that mesic forests can be found at higher elevations, on steeper slopes and there is a lower proportion of bedrock on the surface. Mesic temperate forests appear at altitude of 1100 m (Matevski et al. 2011), but they can appear at lower altitudes in a special geomorphological conditions, as ravines and valleys Čarni et al. | Warm and mesic forests and this lower the explanatory power of the model. Since no precise digital elevation model exists, we used slope and aspect as such, not considering other geomorphological features (Jiménez-Alfaro et al. 2014a; Dirnboeck et al. 2003). The second model comprises ecological variables (Tables 1 and 2; Fig. 3 and supplementary Appendix 1). It is evident that temperatures are higher in the lowlands and there is more nutrient in the soil due to quicker decomposition. A certain caution is needed with this model, since indicator values are calculated on the basis of species composition and often show too high correlation with the data (Zelený and Schaffers 2012). The third model is based on structural and functional characteristics (Tables 1 and 2, Fig. 3, supplementary Appendix 1). Higher structural and functional diversity in the lowlands can be observed: more species, trees, scrubs, more chamaephytes and therophytes. The canopy is closer at higher altitudes in mesic forest. Higher structural and functional diversity of lowland warm temperate forest is also due to the higher human impact (cutting, grazing etc.) that lasted in the area since centuries (Matevski et al. 2008). Their explanatory power is high but below that of the models Ecology and Chorotypes (Pellisser et al. 2010). The fourth model was built on the basis of the content of K, pH, C/N, sand, Mg and clay (Tables 1 and 2, Fig. 3, supplementary Appendix 1). Only pH shows a positive linear correlation with AX1. pH is lower at high altitudes, since decomposition is slower and processes of acidifications due to higher precipitation are more intense; bedrock does not therefore have such an influence as in lowlands, where decomposition is much faster (Filipovski 1995; Zilioli et al. 2011). The fifth model was based on Chorotypes (Tables 1 and 2, Fig. 3, supplementary Appendix 1). The results are in line with the previous models and show that more species originating from the north (Eurasian and boreal species), as well as South European oreophytes can be found in mesic temperate forests (Čarni et al. 2009; Marinšek et al. 2013). The sixth model comprises all variables. As expected this model has the best explanatory power. The best model comprises some variables of all five models that show that all variables contribute to the diversity of vegetation. Comparison of the six models (Fig. 4) shows that the most important factor is ecological conditions. The reason for the high explanatory power of this set of variables has already been discussed (Zelený and Schaffers 2012). Because the development of vegetation since the Pleistocene is one of the main factors that determines the actual vegetation of the Earth (RivasMartínez et al. 2011), the high explanatory power of chorotypes, representing the origin and development of flora and vegetation, was expected. However, it is surprisingly high in comparison to the model Structure, a combination of structural and functional characteristics. The latter are based on life and growth forms, which are among the most well-known plant traits associated with ecosystem functions and revealing plant strategies related to both biotic and abiotic factors (Wellstein et al. 2013). But at the same time, we should also consider the origin of flora and evolution during and after Pleistocene (see further). On the other side, the reason might be also the relatively limited research area, with unique environmental conditions and fairly diverse origin of the flora and vegetation. The low explanatory value of soil characteristics is to be expected, because this set of variables seems to be the most conservative among all sets (Paušič and Čarni 2013). In addition to the origin of flora and the past development of vegetation, horotypes reflect the evolutionary aspect. For instance, Fagus sylvatica, one of the dominant species of mesic forests, once treated as Fagus moesica (Czeczott 1933), shows a certain degree of genetic peculiarity (Brus 2010). Research of this kind of variability is hard to detect and there are hardly any paper dealing with evolution on ecological scale in the region. But we can find some evidences on phenotypic diversification. For instance, European souslik (Spermophilus citellus) a common animal species of the lowlands of Eastern Europe, appears in a phenetically well-differentiated population in high altitudes in the central part of the Balkan Peninsula and may be a relict of the late Pleistocene, when expanding forests restricted it to a few high mountain pastures, where it diverged allopatically (Kryštrufek 1996). We can also detect a diversification within the lineage of Cistus reflected in plant traits due to adaptive radiation (Čarni et al. 2010; Guzmán et al. 2009; Onstein et al. 2014). During last years have appeared an increasing number of research papers about rapid evolutionary response to Downloaded from https://academic.oup.com/jpe/article/9/4/421/2222523 by guest on 03 June 2022 Figure 4: comparison of accuracy of the models (their explanatory power) for AX1. The thick full line represents the function for real data, to which we added result functions for each individual model. The chart clearly indicates the models’ explanatory power, presented in the diagram with a function incline (k) with the formula y = k(x) + n. 429 430 that remained within the region during the Pleistocene (Petit et al. 2008; Rull 2010). These refugia are often too localized or too small to be detected by pollen analysis or other conventional paleoecological techniques (Birks and Willis 2008). The high diversity of chorotypes suggests that macorefugia for warm temperate forests and also cryptic refugia of mesic temperate forests existed in the region. The macrorefugia were more species rich but cryptic refugia also had their floristic individuality (Birks and Willis 2008). Both are reflected in the present vegetation of the region. In the rather small region and under stable and fairly uniform ecological conditions, chorotypes better explain the transition between warm and mesic temperate forests than the model based on structure and function. CoNClusIoN The research took place in the Balkan Mountains where the Pleistocene climatic oscillations were not so pronounced as in other parts of the Europe. The results show that the model based on all variables has the best explanatory power and it is followed by models based on groups of variables, as: ecology, chorotypes, geomorphology, structure and function and soil properties. The relatively high explanatory power of chorotypes supports the idea about eco-evolutionary dynamics on short ecological time-scale. In the future the results should be tested also in regions with major climatic dynamics. suPPlEmENtarY matErIal Supplementary material is available at Journal of Plant Ecology online. FuNDINg National Park Galičica, Research Centre of the Slovenian Academy of Sciences and Arts, Slovenian Research Agency (P1-0236). aCKNoWlEDgEmENts For help during the field work, we are grateful to director Zoran Angeloski and the staff of the Public Institution Galičica National Park and Til Dieterich (GFA Consulting Group), who provided us with logistical support during the field work. The soil analyses were performed by the Pedological Laboratory at the Biotechnical Faculty, University of Ljubljana. We are grateful to Andrej Blejec for useful advice and comments during statistical elaboration of the data. Thanks to Iztok Sajko for prepering the figures for printing. We also thank two anonymous reviewers for their useful comments and suggestions. Conflict of interest statement. None declared. rEFErENCEs Acevski J (2001) Dendroflorata na planinata Galičica. Doktorska disertacija: Šumarski fakultet. Downloaded from https://academic.oup.com/jpe/article/9/4/421/2222523 by guest on 03 June 2022 environmental change (Broennimann et al. 2007; Buswell et al. 2011). We can thus hypothesise the possibility of the beginning of allopatric speciation in the Balkans, although further genetic analysis should be done to confirm such an assumption. In the past, relationship between evolutionary and ecological processes were decoupling, since the evolutionary processes were understood as acting on a much longer scale that ecological ones shaping communities and vegetation. But in last decades a growing realization emerged that evolution occurs on ecological time scales and this should be accounted for ecological research (Carroll et al. 2007; Eriksson 2013, 2014; Qian and Jiang 2014). Changing of climate after Pleistocene with new ecological circumstances could provoke an ecological niche shift of species or even additional evolutionary response affecting the fundamental niche (Eriksson 2013; Pearman et al. 2008). In this regard also, the species and trait distribution cannot be considered as constant but as dynamic and changing with regard to their niche relationships (Pearman et al. 2008). If we compare mesic temperate forest in the broader region of southeast-Europe, we can find much larger differences in their structural characteristics evaluated by traits then in their floristic composition (Marinšek et al. 2013; Tzonev et al. 2006). In this way, the species in the region might adapt to the new ecological circumstances, maintaining their origin, but had changed their niche and also traits consequently. The development of vegetation since the Pleistocene is crucial in the establishment of the region’s species pool that builds the todays vegetation (Zelený 2009). The main refugia of beech forests can be found in the north-western Balkans (Magri et al. 2006; Médail and Diadema 2009; Marinšek et al. 2013; Willner et al. 2009), in the region beyond study area. On the other hand, refugia of warm temperate forests can be found in southern parts of the Balkans (Petit et al. 2002). It seems that drought resistant species survived in southEuropean glacial refugia, whereas more mesophilous ones disappeared (Michelot et al. 2012; Svenning 2003). Petit et al. (2005) suggest that most tree species of mesic temperate Euro-Asian flora survived in the north of the Mediterranean region (higher than 45° latitude). In our case this would mean north-western outcrops of the Dinaric Alps (Brus 2010). Such anticipating is supported also by floristic analysis of mesophilous beech forests of the whole Balkan Peninsula. The flora of beech forests in the north-western Balkans is much richer than in the central and southern Balkans. At the same time, both groups have their own diagnostic species (Marinšek et al. 2013). As the floristic inventory of individual sample plots (relevés) in higher in warm temperate forests, where more than 41 plant species on average can be found on a single plot, whereas only 25 can be found in mesic temperate forests (supplementary Appendix 1). All these facts allow us to look for the cradle of warm temperate forest in the region (Petit et al. 2002). However, the high endemism in the region, especially of mesic temperate forests, and also genetic variability can be explained only by cryptic refugia of mesophilous vegetation Journal of Plant Ecology Čarni et al. | Warm and mesic forests 431 Adams J (2009) Species Richness. Patterns in Diversity of Life. Chichester, UK: Praxis publishing. Aertsen W, Kint V, van Orshoven J, et al. (2010) Comparison and ranking of different modelling techniques for prediction of site index in Mediterranean mountain forests. Ecol Model 221:1119–30. Akaike H (1973) Information theory and an extension of the maximum likelihood principle. In Petran BN, Csari F (eds). International Symposium on Information Theory, 2nd edn. Budapest: Akademiai Kiado. Alahuhta J, Heino J, Luoto M (2011) Climate change and the future distributions of aquatic macrophytes across boreal catchments. J Biogeogr 38:383–93. Cubizolle H, Argant J, Fassion F, et al. (2014) Vegetation history from the end of the late-glacial and human impact from the midHolocene in the Eastern Massif Central (France). Quaternaire 25:209–36. Czeczott H (1933) A study on the variability of the leaves of beeches: F. orientalis Lipsky, F. sylvatica L. and intermediate forms. Part I. Rocznik Dendrologiczny 5:45–121. Danz NP, Frelich LE, Reich PB, et al. (2013) Do vegetation boundaries display smooth or abrupt spatial transitions along environmental gradients? Evidence from the prairie-forest biome boundary of historic Minnesota, USA. J Veg Sci 24:1129–40. Aronof S (2005) Remote Sensing for GIS Managers. New York, NY: ESRI Press. Diklić N (1986) Flora SR Srbije, Vol. X. Beograd: Srpska akademija nauka i umetnosti. Bergmeier E, Dimopoulos P (2008) Identifying plant communities of thermophilous deciduous forest in Greece: species composition, distribution, ecology and syntaxonomy. Plant Biosyst 142:228–54. Dirnboeck T, Dullinger S, Gottfried M, et al. (2003) Mapping alpine vegetation based on image analysis, topographic variables and Canonical Correspondence Analysis. Appl Veg Sci 6:85–96. Birks HJB, Willis KJ (2008) Alpines, trees, and refugia in Europe. Plant Ecol Divers 1:147–60. Ewald J (2003a) The calcareous riddle: why are there so many calciphilous species in the Central European flora? Folia Geobotanica 38:357–66. Braun-Blanquet J (1964) Pflanzensoziologie. Vegetationskunde. Wien: Springer. Grundzüge der Broennimann O, Treier UA, Müller-Schärer H, et al. (2007) Evidence of climatic niche shift during biological invasion. Ecology Lett 10:701–9. Brus R (2010) Growing evidence for the existence of glacial refugia of European beech (Fagus sylvatica L.) in the south-eastern Alps and north-western Dinaric Alps. Periodicum Biologorum 112:239–46. Buswell JM, Moles AT, Hartley S (2011) Is rapid evolution common in introduced plant species? J Ecol 99:214–24. Burnham KP, Anderson DR (2002) Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. Vienna: Springer. Campetella G, Botta-Dukát Z, Wellstein C, et al. (2011) Patterns of plant trait–environment relationships along a forest succession chronosequence. Agric Ecosyst Environ 145:38–48. Caprio E, Ellena I, Rolando A (2009) Assessing habitat/landscape predictors of bird diversity in managed deciduous forests: a seasonal and guild-based approach. Biodivers Conserv 18:1287–303. Carroll SP, Hendry AP, Reznick DN, et al. (2007) Evolution on ecological time‐scales. Funct Ecol 21:387–93. Čarni A, Matevski V, Šilc U (2010) Morphological, chorological and ecological plasticity of Cistus incanus in the southern Balkans. Plant Biosyst 144:602–17. Čarni A, Košir P, Karadžić B, et al. (2009) Thermophilous deciduous forests in southeastern Europe. Plant Biosyst 143:1–13. Čarni A, Matevski V (2015) Impact of climate change on mountain flora and vegetation in the Republic of Macedonia (Central part of the Balkan Peninsula). In Ozturk M, Hakeem KR, Faridah-Hanum I, Efe R (eds). Climate Change Impacts on High-Altitude Ecosystems. Heidelberg: Springer, 189–214. Chytrý M, Tichý L, Roleček J (2003) Local and regional patterns of species richness in Central European vegetation types along the pH/calcium gradient. Folia Geobotanica 38:429–42. Cornelissen JHC, Lavorel S, Garnier E, et al. (2003) A handbook of protocols for standardised and easy measurement of plant functional traits worldwide. Aust J Bot 51:335–80. Ewald J (2003b) The sensivity of Ellenberg indicator values to the completeness of vegetation relevés. Basic Appl Ecol 4:507–13. Eriksson O (2013) Species pools in cultural landscapes–niche construction, ecological opportunity and niche shifts. Ecography 36:403–13. Eriksson O (2014) Vegetation change and eco-evolutionary dinamics. J Veg Sci 25:1141–7. Faraway JJ (2005) Linear Models with R. Boca Raton, FL: Chapman & Hall/CRC. Feliner GF (2011) Southern European glacial refugia: a tale of tales. Taxon 60:365–72. Feliner GN (2014) Patterns and processes in plant phylogeography in the Mediterranean Basin. A review. Persp Plant Ecol Evol Syst 16:265–78. Filipovski G (1995) Počvite na Republika Makedonija. Skopje: MANU. Filipovski G, Rizovski P, Ristevski P (1994) Karakteristiki na klimatskovegetacisko-počvenite zoni (regioni) vo Republika Makedonija. Skopje: MANU. Frajman B, Oxelman B (2007) Reticulate phylogenetics and phytogeographical structure of Heliosperma (Sileneae, Caryophyllaceae) inferred from chloroplast and nuclear DNA sequences. Mol Phylogenet Evol 43:140–55. Gajić M (1980) Pregled vrsta flore SR Srbije sa biljnogeografskim oznakama. Glasnik Šumarskog fakulteta, Beograd A 54:111–41. Guzmán B, Lledó MD, Vargas P (2009) Adaptive radiation in mediterranean Cistus (Cistaceae). PLOS ONE 4:e6362. Halder A, Ghosh A, Ghosh S (2011) Supervised and unsupervised landuse map generation from remotely sensed images using ant based systems. Appl Soft Comput 11:5770–81. Haslem A, Kelly LT, Nimmo DG, et al. (2011) Habitat or fuel? Implications of long-term, post-fire dynamics for the development of key resources for fauna and fire. J Appl Ecol 48:247–56. Hegi G (1946–2003) Illustrierte Flora von Mitteleuropa, Vol. I–VI. Berlin: Paul Parey Verlag. Downloaded from https://academic.oup.com/jpe/article/9/4/421/2222523 by guest on 03 June 2022 Armas C, Rodríguez-Echeverría S, Pugnaire FI (2011) A field test of the stress-gradient hypothesis along an aridity gradient. J Veg Sci 22:1–10. Dawson CW, Abrahart RJ, See LM (2007) Hydrotest: a web-based toolbox of evaluation metrics for the standardised assessment of hydrological forecasts. Environ Model Softw 22:1034–52. 432 Journal of Plant Ecology Hennekens SM, Schaminée JHJ (2001) TURBOVEG, a comprehensive database management systeme for vegetation data. J Veg Sci 12:589–91. Matevski V, Čarni A, Avramoski O, et al. (2011) Forest Vegetation of the Galičica Mountain Range in Macedonia. Ljubljana: ZRC-Publishing & JU Nacionalen park Galičica. ISO 10390 (1996) Soil quality – Determination of pH:5. Médail F, Diadema K (2009) Glacial refugia influence plant diversity pattern in the Mediterranean basin. J Biogeogr 36:1333–45. ISO 11261 (1996) Soil quality – Determination of total nitrogen – Modified Kjeldahlmethod:4. ISO 11464 (1996) Soil quality – Pretreatment of samples for physicochemical analyses:9. ISO 14235 (1999) Soil quality – Determination of organic carbon by sulphocromic oxidation:5. Jiménez-Alfaro B, Marceno C, Bueno A, et al. (2014a) Biogeographic deconstruction of alpine plant communities along altitudinal and topographic gradient. J Veg Sci 25:160–71. Jordanov D (1963–1982) Flora na Narodna Republika Bulgarija, Vol. I–VIII. Sofia: Bulgarskata akademija na naukite. Josifović M (1970–1977) Flora SR Srbije, Vol. I–IX. Beograd: Srpska akademija nauka i umetnosti. Košir P, Čarni A, Di Pietro R (2008) Classification and phytogeographical differentiation of broad-leaved ravine forests in southeastern Europe. J Veg Sci 19:331–42. Michelot A, Bréda N, Damesin C, et al. (2012) Differing growth responses to climatic variations and soil water deficiency of Fagus sylvatica, Quercus petraea and Pinus sylvestris in a temperate forests. For Ecol Manag 265:161–71. Nie L, Chen Y, Chu H (2011) Asymptotic variances of maximum likelihood estimator for the correlation coefficient from a BVN distribution with one variable subject to censoring. J Stat Plan Infer 141:392–401. Oksanen J, Kindt R, Legendre P, et al. (2013) Vegan-package. http:// cran.r-project.org, http://vegan.r-forge.r-project.org/ (27 March 2015, date last accessed). Onstein RE, Carter RJ, Xing Y, et al. (2014) Diversification rate shifts in the Cape Floristic Region: the right traits in the right place at the right time. Persp Plant Ecol Evol Syst 16:331–40. Pearman PB, Guisan A, Broennimann O, et al. (2008) Niche dynamics in space and time. Trends Ecol Evol 23:149–58. Košir P, Cassavechia S, Čarni A, et al. (2013) Ecological and phytogeographical differentiation of oak-hornbeam forests in southeastern Europe. Plant Biosyst 147:84–98. Pearson K (1896) Mathematical contributions to the theory of evolution. III. Regresion, heredity and panmixia. Philos Trans Roy Soc Lond 187:253–318. Kraft N J, Godoy O, Levine J M (2015) Plant functional traits and the multidimensional nature of species coexistence. Proc Natl Acad Sci USA 112:797–802. Pellisser L, Fournier B, Guisan A (2010) Plant traits co-vary with altitude in grasslands and forests in the European Alps. Plant Ecol 211:351–65. Kropf M, Comes HP, Kadereit JW (2008) Causes of the genetic architecture of south-west European high mountain disjuncts. Plant Ecol Divers 1:217–28. Petit RJ, Brewer S, Bordács S, et al. (2002) Identification of refugia and post-glacial colonisation routes of European white oaks based on chloroplast DNA and fossil pollen evidence. For Ecol Manag 156:49–74. Kryštrufek B (1996) Phenetic variation in the European souslik, Spermophilus citellus (Mammalia: Rodentia). Bonner Zoologische Beiträge 46:93–109. Kučera J, Marhold K, Lihová J (2010) Cardamine maritima group (Brassicaceae) in the amphi-Adriatic area: a hotspot of species diversity revealed by DNA sequences and morphological variation. Taxon 59:148–64. Petit RJ, Hampe A, Cheddadi R (2005) Climate changes and tree phylogeography in the Mediterranean. Taxon 54:877–85. Petit RJ, Hu FS, Dick CW (2008) Forest of the past: a window to future changes. Science 320:450–2. Pignatti S (1982) Flora d’Italia, Vol. 1. Bologna: Edagricole. Lazarevski A (1993) Klimata vo Makedonija. Skopje: Kultura. Pignatti S, Menegoni P, Pietrosanti S (2005) Valori di bioindicazione delle piante vascolari della flora d`Italia. Braun-Blanquetia 39:1–97. Lang P, Ewald J (2014) Predictive modelling and monitoring of Ellenberg moisture value validates restoration success in floodplain forests. Appl Veg Sci 17:543–55. Qian H, Jiang L (2014) Phylogenetic community ecology: integrating community ecology and evolutionary biology. J Plant Ecol 7:97–100. Lehmann A, Overton J, Leathwick J (2002) GRASP: generalized regression analysis and spatial prediction. Ecol Model 157:189–207. Quézel P, Médail F (2003) Ecologie et Biogéographie des Forêts du Bassin Méditerranéen, Vol. 572. Paris, France: Elsevier. Magri D, Vendramin GG, Comps B, et al. (2006) A new scenario for the Quaternary history of European beech populations: palaeobotanical evidence and genetic consequences. New Phytol 171:199–221. Raunkiaer C (1934) The Life Forms of Plants and Statistical Plant Geography. Oxford, UK: Charendon Press. Marinšek A, Šilc U, Čarni A (2013) Geographical and ecological differentiation of mesophilous Fagus forest vegetation in the Southeast Europe. Appl Veg Sci 13:131–47. Matevski V, Čarni A, Kostadinovski M, et al. (2010) Notes on pyhtosociology of Juniperus excelsa in Macedonia (Southern Balkan Peninsula). Hacquetia 9:161–5. Matevski V, Čarni A, Kostadinovski M, et al. (2008) Flora and Vegetation of the Macedonian Steppe. ZRC – Publishing. Paušič A, Čarni A (2013) Records of past land use are best stored in soil properties. Plant Biosyst 147:654–63. Reger B, Kölling C, Ewald J (2011) Modeling effective thermal climate for mountain forests in the Bavarian Alps: which is the best model? J Veg Sci 22:677–87. Ricklefs RE (1987) Community diversity: relative roles of local and regional processes. Science 235:167–71. Rivas-Martínez S, Rivas Sáenz S, Penas A (2011) Worldwide bioclimatic classification system. Glob Bot 1:1–634. Downloaded from https://academic.oup.com/jpe/article/9/4/421/2222523 by guest on 03 June 2022 Jiménez-Alfaro B, Chytrý M, Rejmánek M, et al. (2014b) The number of vegetation types in European countries: major determinants and extrapolation to other regions. J Veg Sci 25:863–72. Micevski K (1985–2005) Flora na Republika Makedonija, Vol. I–IV. Skopje: MANU. Čarni et al. | Warm and mesic forests Ruiz-Labourdette D, Nogués-Bravo D, Sáinz Ollero H, et al. (2012) Forest composition in Mediterranean mountains is projected to shift along the entire elevation gradient under climate change. J Biogeogr 39:162–76. Rull V (2010) On microrefugia and cryptic refugia. J Biogeogr 37:1623–5. 433 Ter Braak CJF, Šmilauer P (2002) CANOCO Reference Manual and CanoDraw for Windows User’s Guide. Software for Canonical Community Ordination (Version 4.5). Ithaca: Microcomputer Power. Tichý L (2002) JUICE, software for vegetation classification. J Veg Sci 13:451–3. Tzonev R, Dimitrov M, Chytrý M, et al. (2006) Beech forest communities in Bulgaria. Phytocoenologia 36:247–79. Wehn S, Lundemo S, Holten JI (2014) Alpine vegetation along multiple environmental gradients and possible consequences of climate change. Alpine Bot 124:155–64. Saatkamp A, Römermann C, Dutoit T (2010) Plant functional traits show non-linear response to grazing. Folia Geobotanica 45:239–52. Wellstein C, Chelli S, Campetella G, et al. (2013) Intraspecific phenotypic variability of plant functional traits in contrasting mountain grasslands habitats. Biodivers Conserv 22:2353–74. Seddon AW, Macias-Fauria M, Willis KJ (2015) Climate and abrupt vegetation change in Northern Europe since the last deglaciation. The Holocene 25:25–36. Willner W, Di Pietro R, Bergmeier E (2009) Phytogeographical evidence for post-glacial dispersal limitation of European beech forest species. Ecography 32:1011–8. StatSoft, Inc. (2011) Electronic Statistics Textbook. Tulsa, US: StatSoft. http://www.statsoft.com/textbook/. Wood S, Augustin N (2002) GAMs with integrated model selection using penalized regression splines and applications to environmental modelling. Ecol Model 157:157–77. Strid A, Tan K (1997) Flora Hellenica, Vol. I. Königstein: `Koeltz Scientific Books. Strid A, Tan K (2002) Flora Hellenica, Vol. II. Ruggell: Gantner Verlag. Surina, B, Schönswetter P, Schneeweiss GM (2011) Quaternary range dynamics of ecologically divergent species (Edraianthus serpyllifolius and E. tenuifolius, Campanulaceae) within the Balkan refugium. J Biogeogr 38:1381–93. Svenning JC (2003) Deterministic Plio-Pleistocene extinctions in the European cool-temperate tree flora. Ecol Lett 6:646–53. Swets K (1988) Measuring the accurancy of diagnostic systems. Science 240:1285–93. Tamstorf MP, Illeris L, Hansen BU, et al. (2007) Spectral measures and mixed models as valuable tools for investigating controls on land surface phenology in high arctic Greenland. BMC Ecol 7:9. Wood SN (2006) Generalized Models. An Introduction with R. Boca Raton, FL: Chapman & Hall/CRC. Zelený D (2009) Co‐occurrence based assessment of species habitat specialization is affected by the size of species pool: reply to Fridley et al.(2007). J Ecol 97:10–7. Zelený D, Schaffers AP (2012) Too good to be truth: pitfalls of using mean Ellenberg indicator values in vegetation analyse. J Veg Sci 23:419–31. Zilioli DM, Bini C, Wahsha M, et al. (2011) The pedological heritage of the Dolomites (Northern Italy): features, distibution and evolution of the soils, with some implications for land management. Geomorphology 135:232–47. Zuur AF, Ieno EN, Walker NJ, et al. (2009) Mixed Effects Models and Extensions in Ecology with R. New York, NY: Springer. Downloaded from https://academic.oup.com/jpe/article/9/4/421/2222523 by guest on 03 June 2022 R Development Core Team (2013) R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. http://www.R-project.org (27 March 2015, date last accessed).