Josep Padullés Cubino

Borja Jiménez-Alfaro

Francesco Maria Sabatini

Wolfgang Willner

Zdeňka Lososová

Idoia Biurrun

Jörg Brunet

Juan Antonio Campos

Adrian Indreica

Florian Jansen

Jonathan Lenoir

Željko Škvorc

Kiril Vassilev

Milan Chytrý

Josep Padullés Cubino

Josep Padullés
email: padullesj@gmail.com

Plant taxonomic and phylogenetic turnover increases toward climatic extremes and depends on historical factors in European beech forests

Keywords: beta-diversity, climatic gradient, community phylogenetics, deciduous forests, distribution range margin, European Vegetation Archive (EVA), Fagus sylvatica, plant diversity, post-glacial dispersal, species rarity, vegetation plots

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

A key issue in macroecology is understanding the role of biogeographical processes in determining spatial variation in biodiversity [START_REF] Ricklefs | A comprehensive framework for global patterns in biodiversity[END_REF]. Classical biogeographical studies have mainly asked why some sites have more species than others, but much less attention has been paid to explaining spatial variation in the composition of communities, i.e., beta-diversity [START_REF] Qian | Beta diversity in relation to dispersal ability for vascular plants in North America[END_REF][START_REF] Kraft | Disentangling the drivers of diversity along latitudinal and elevational gradients[END_REF][START_REF] Myers | Beta-diversity in temperate and tropical forests reflects dissimilar mechanisms of community assembly[END_REF]. Exploring how betadiversity varies among regions with contrasting ecological and evolutionary histories is fundamental for understanding the processes that structure natural communities [START_REF] Leibold | The metacommunity concept: a framework for multi-scale community ecology[END_REF] and for guiding conservation efforts [START_REF] Socolar | How should beta-diversity inform biodiversity conservation[END_REF].

Environmental, spatial, and historical factors can all drive patterns in beta-diversity.

First, environmental filtering and species' ecological requirements (niche) determine species' presence in a given habitat [START_REF] Macarthur | Geographical ecology[END_REF][START_REF] Tilman | Plant strategies and the dynamics and structure of plant communities[END_REF]. Sites with more dissimilar environments should host more different sets of species and thus exhibit greater betadiversity [START_REF] Whittaker | Communities and ecosystems[END_REF]. Second, spatial processes can also influence species composition through stochastic events, ecological drift, and species' distinct dispersal abilities [START_REF] Hubbell | The unified neutral theory of biodiversity and biogeography[END_REF][START_REF] Nekola | The distance decay of similarity in biogeography and ecology[END_REF][START_REF] Leibold | The metacommunity concept: a framework for multi-scale community ecology[END_REF][START_REF] Soininen | The distance decay of similarity in ecological communities[END_REF]. Third, historical processes, including speciation, extinction, range expansion from areas of origin, or recolonization from refugia, are also important drivers of species distributions (Svenning et al., 2015). Yet, it remains poorly understood how these three sets of drivers combine to shape beta-diversity across large spatial scales, and in specific habitat types.

The relationships between environmental, spatial, and historical processes and betadiversity have often been examined using taxonomic information only (e.g., [START_REF] Qian | Beta diversity in relation to dispersal ability for vascular plants in North America[END_REF][START_REF] Keil | Patterns of beta diversity in Europe: the role of climate, land cover and distance across scales[END_REF][START_REF] Sabatini | Betadiversity of central European forests decreases along an elevational gradient due to the variation in local community assembly processes[END_REF]. However, species are not evolutionarily independent, and historical processes coupled with species diversification are expected to leave significant imprints on communities' phylogenetic structure [START_REF] Emerson | Phylogenetic analysis of community assembly and structure over space and time[END_REF]. Combining taxonomic and phylogenetic measures of beta-diversity can help understand the effect of ecological and evolutionary processes on biodiversity patterns [START_REF] Graham | Phylogenetic beta diversity: linking ecological and evolutionary processes across space in time[END_REF]Swenson, 2011a). For example, in areas with a large fraction of narrow-ranging species with recent divergence time (i.e., neo-endemics), taxonomic beta-diversity should be high but phylogenetic beta-diversity low (Swenson, 2011a). In contrast, if the same narrow-ranging species have relatively old divergence time (i.e., belong to lineages with long-standing and disparate evolutionary histories), both taxonomic and phylogenetic beta-diversity should be high [START_REF] Graham | Phylogenetic beta diversity: linking ecological and evolutionary processes across space in time[END_REF]). Yet, taxonomic and phylogenetic beta-diversities have seldom been employed over large spatial extents to disentangle the relative importance of ecological and historical processes (Swenson, 2011b;[START_REF] Hardy | Phylogenetic turnover in tropical tree communities: impact of environmental filtering, biogeography and mesoclimatic niche conservatism[END_REF][START_REF] Kubota | Beta diversity of woody plants in the Japanese archipelago: the roles of geohistorical and ecological processes[END_REF].

Several features make European beech forests a good model system to test the influence of environmental and historical factors on the patterns of taxonomic and phylogenetic beta-diversity. European beech (Fagus sylvatica) is one of the most abundant deciduous tree species in Europe, which typically forms monospecific, dense canopies, with a relatively homogeneous understory. It extends from the mountains of southern Europe (northern Spain, Italy, and Greece) to the lowlands in south-eastern England and southern Sweden [START_REF] Bolte | The north-eastern distribution range of European beecha review[END_REF][START_REF] Caudullo | Chorological maps for the main European woody species[END_REF]. Its distribution is mainly constrained by dry summers and frosty winters [START_REF] Houston Durrant | Fagus sylvatica in Europe: distribution, habitat, usage and threats[END_REF]. Nonetheless, many of its understory species have different ecological and geographical requirements than beech itself and may be restricted to specific regions and ecological conditions [START_REF] Willner | Phytogeographical evidence for postglacial dispersal limitation of European beech forest species[END_REF]. Furthermore, the geographical distribution of beech forest species has repeatedly been shown to be affected by dispersal limitations partially reflecting the history of post-glacial recolonization from refugia [START_REF] Magri | Patterns of post-glacial spread and the extent of glacial refugia of European beech (Fagus sylvatica)[END_REF][START_REF] Lehsten | Modelling the Holocene migrational dynamics of Fagus sylvatica L. and Picea abies (L.) H. Karst: LPJ-CATS migration model[END_REF][START_REF] Willner | Phytogeographical evidence for postglacial dispersal limitation of European beech forest species[END_REF]. Several studies have examined beech forest diversity [START_REF] Willner | Classification of European beech forests: a Gordian Knot?[END_REF] and its environmental and historical drivers at the regional [START_REF] Hrivnák | Species richness pattern along altitudinal gradient in central European beech forests[END_REF][START_REF] Ujházyová | Diversity of beech forest vegetation in the Eastern Alps, Bohemian Massif and the Western Carpathians[END_REF][START_REF] Weigel | Forest understory vegetation is more related to soil than to climate towards the cold distribution margin of European beech[END_REF] and continental scale [START_REF] Jiménez-Alfaro | History and environment shape species pools and community diversity in European beech forests[END_REF][START_REF] Willner | Alpha and beta diversity in central European beech forests[END_REF][START_REF] Willner | Phytogeographical evidence for postglacial dispersal limitation of European beech forest species[END_REF].

However, no studies have simultaneously tested for the effect of these drivers on both taxonomic and phylogenetic turnover (i.e., the beta-diversity component that reflects the change in species identities and phylogenetic relationships between two or more assemblages; [START_REF] Baselga | Partitioning the turnover and nestedness components of beta diversity[END_REF] across the whole range of beech forests.

Here, we used plant occurrences from vegetation-plot records and aggregated them within grid cells of 0.25º × 0.25º spatial resolution across the range of European beech forests. Then we calculated the spatial turnover component of taxonomic and phylogenetic beta-diversity between a given cell and all the other cells. Specifically, we tested two hypotheses. First, (H1) taxonomic and phylogenetic spatial turnover increases in the cells with more distinct environmental conditions or closer to the margin of the geographical distribution of European beech [START_REF] Macarthur | Geographical ecology[END_REF][START_REF] Whittaker | Communities and ecosystems[END_REF]. We expected this because beech forests closer to the edges of their environmental and geographical distribution can suffer from stronger climatic and environmental stress, and host unique sets of species, some of them potentially deriving from other habitat types. Second, (H2) we expected phylogenetic turnover to be low in relation to taxonomic turnover when moving closer to the putative glacial refugial areas of beech. These areas, mostly located on three southern European peninsulas (Iberian, Italian and Balkan) [START_REF] Magri | A new scenario for the Quaternary history of European beech populations: palaeobotanical evidence and genetic consequences[END_REF], can host a higher fraction of narrow-ranging and rare species within the understory of beech forests. In some cases, such species are closely related because they evolved as a consequence of geographical isolation during the cold stage (i.e., by vicariance) and did not disperse into newly available environments after glaciation [START_REF] Nekola | Paleorefugia and neorefugia: the influence of colonization history on community patterns and process[END_REF][START_REF] Turner | Explaining the global biodiversity gradient: energy, area, history and natural selection[END_REF][START_REF] Svenning | Could the tree diversity pattern in Europe be generated by postglacial dispersal limitation[END_REF].

Methods

Vegetation data

We compiled georeferenced vegetation-plot records from the European Vegetation Archive (EVA; [START_REF] Chytrý | European Vegetation Archive (EVA): an integrated database of European vegetation plots[END_REF]accessed on June 28th, 2019). This database contains more than 1.5 million vegetation plots sampled across Europe (see Appendix S1 for an overview of contributing databases). European beech (Fagus sylvatica) forests are amongst the bestsampled vegetation types in EVA (~ 35,000 vegetation plots), providing a unique model to address macroecological questions [START_REF] Jiménez-Alfaro | History and environment shape species pools and community diversity in European beech forests[END_REF]. To select vegetation plots dominated by European beech forest from the EVA database, we used a classification expert system of European EUNIS habitats [START_REF] Chytrý | EUNIS Habitat Classification: Expert system, characteristic species combinations and distribution maps of European habitats[END_REF] run in the JUICE program [START_REF] Tichý | JUICE, software for vegetation classification[END_REF]. This expert system contains definitions of individual EUNIS habitats based on their species composition and geographic location [START_REF] Chytrý | EUNIS Habitat Classification: Expert system, characteristic species combinations and distribution maps of European habitats[END_REF]. We focused only on European beech forests to reduce the potential confounding effects of different evolutionary histories and species pools across habitat types.

We selected plots sampled after 1970 with areas ranging between 100 and 1,000 m 2 , which were the most common plot sizes in the database, and avoided plots in which the location uncertainty of coordinates was larger than 10 km. However, in specific regions with lower sampling effort, we retained plots with unknown location uncertainty. For each plot, we obtained data on species occurrences and selected only angiosperms, thus excluding gymnosperms, pteridophytes, lycophytes, and non-vascular plants (~4% of total species). We focused on angiosperms to avoid inflated results due to a few taxa connected to deep phylogenetic nodes [START_REF] Massante | Environment and evolutionary history depict phylogenetic alpha and beta diversity in the Atlantic coastal white-sand woodlands[END_REF][START_REF] Qian | Geographic patterns and environmental correlates of taxonomic and phylogenetic beta diversity for large-scale angiosperm assemblages in China[END_REF] see Appendix S2 for further justification and results including all vascular plants). We standardized the species names in our dataset according to The Plant List using the R package 'Taxonstand' [START_REF] Cayuela | Taxonstand: Taxonomic Standardization of Plant Species Names. R package version 2[END_REF].

We aggregated vegetation plots in a raster of grid cells of 0.25º × 0.25º spatial resolution and derived the composition of the species pool of each cell (N = 2,035; total number of species = 2,213). The species pool of each cell was obtained by aggregating all plant occurrences, independent of the vegetation layer (understory and upper canopy layers), from the final set of vegetation plots classified as beech habitat and contained within the focal cell. Grid cells with less than five vegetation plots (1,049) were removed from subsequent analyses to guarantee a reliable identification of the pool of species found in the cell (N = 986; total number of species = 2,111). To avoid effects of a disproportionally high density of sample plots in particular grid cells, we conducted a heterogeneity-constrained random (HCR) resampling [START_REF] Lengyel | Heterogeneity-constrained random resampling of phytosociological databases[END_REF] with the R package 'vegclust' [START_REF] De Cáceres | The management of vegetation classifications with fuzzy clustering[END_REF]. We did this separately for each grid cell and selected the same number (N = 5) of plots per grid cell. The HCR resampling strategy maximizes the mean and minimizes the variance of the compositional dissimilarity between pairs of plots. Consequently, it tends to select plots that are representative of compositional diversity of vegetation within each grid cell. We calculated dissimilarity between vegetation plots co-occurring in the same grid cell with the Simpson dissimilarity index [START_REF] Baselga | Partitioning the turnover and nestedness components of beta diversity[END_REF], computing 1,000 iterations within each grid cell. We also repeated the HCR resampling with 10 and 15 plots per grid cell, but these variants all yielded very similar results to those based on the threshold of five plots per grid cell (Appendix S3).

The total number of vegetation plots included in our study after the HCR resampling was 4,930. They were distributed across 986 grid cells and hosted a total of 1,739 angiosperm species.

Phylogenetic data

We generated a phylogeny by linking the species in our dataset (n = 1,739) with those in the mega-phylogeny implemented in the R package 'V.PhyloMaker' (Jin and [START_REF] Qian | PhyloMaker: an R package that can generate very large phylogenies for vascular plants[END_REF], which includes 74,533 species derived from two mega-trees [START_REF] Zanne | Three keys to the radiation of angiosperms into freezing environments[END_REF][START_REF] Smith | Constructing a broadly inclusive seed plant phylogeny[END_REF]. We used the 'Scenario 3' approach implemented in the same R package to add missing species (~ 27%) to the phylogeny. In this scenario, the tip for a new genus is bound to 1/2 of the family branch, and the new tip of an existing genus is bound to the basal node of a given genus (for further details see [START_REF] Qian | An updated megaphylogeny of plants, a tool for generating plant phylogenies and an analysis of phylogenetic community structure[END_REF][START_REF] Qian | PhyloMaker: an R package that can generate very large phylogenies for vascular plants[END_REF].

Environmental and historical data

We obtained current climatic data from the full set of 19 bioclimatic variables taken from the WorldClim v.2.1 database [START_REF] Fick | WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas[END_REF] at a resolution of 30 seconds (~1 km 2 ) (Appendix S4). We extracted climatic data for each cell with the R package 'raster' [START_REF] Hijmans | raster: Geographic Data Analysis and Modeling[END_REF] and calculated the mean value of all raster pixels contained in the cell. We then used principal component analysis (PCA) to reduce the dimensionality and collinearity of the 19 bioclimatic variables. We standardized all variables to a mean of 0 and a standard deviation of 1 before implementing the PCA. We interpreted the first PCA component (37% of variation explained) as a gradient of aridity, as it distinguished between cells with low temperatures and high precipitation (mainly mountain areas of Central Europe) and cells with high temperatures and low precipitation (mainly Southern Europe; Appendix S4: Table S4.1 and Fig. S4.1). We only used this PCA axis in subsequent analyses for two reasons: (1) aridity has long been recognized as a major factor determining the distribution range of beech forests in Europe [START_REF] Houston Durrant | Fagus sylvatica in Europe: distribution, habitat, usage and threats[END_REF], and ( 2) it showed the strongest relationship with our set of response variables (Appendix S4: Fig. S4.2). We thus refer to the first PCA component as an 'aridity' factor throughout the manuscript.

We additionally downloaded data for soil pH (at 15 cm depth) at a spatial resolution of ~250 m from the SoilGrids database (https://soilgrids.org/; [START_REF] Hengl | SoilGrids250m: Global gridded soil information based on machine learning[END_REF] and calculated the mean value of all raster pixels contained in each cell. We also obtained a shapefile with the geographical distribution of European beech (Fagus sylvatica) from [START_REF] Caudullo | Chorological maps for the main European woody species[END_REF], and used it to calculate the distance (i.e., great-circle distance) between each grid cell and the closest margin of the geographical distribution of European beech (Appendix S5: Fig. S5.1). We calculated distances with the R package 'geosphere' [START_REF] Hijmans | raster: Geographic Data Analysis and Modeling[END_REF]. We assigned zero values to cells located outside the distribution range (~ 9%). Grid cells located closer to the margin of the geographical distribution of European beech were associated with more extreme climatic conditions (Appendix S6: Fig. S6.2), supporting our assumption that geographically marginal areas can suffer from stronger climatic and environmental stress and present greater species turnover. For this reason, we refer to the distance from the distribution margin of beech as an 'environmental variable' througouth the study.

Finally, we also obtained tentative locations of the putative glacial refugia of beech from [START_REF] Magri | A new scenario for the Quaternary history of European beech populations: palaeobotanical evidence and genetic consequences[END_REF]

Species richness (SR) and Species rarity index (SRar)

Taxonomic and phylogenetic beta-diversity might be affected by both the number of species in a cell and their degree of uniqueness or rarity across grid cells. Therefore, we calculated both the total species richness (SR), and a 'species rarity index' (SRar) for each focal grid cell, to explicitly control for these effects. SRar was calculated as follows:

𝑆𝑅𝑎𝑟 = 1 -𝑚𝑒𝑑𝑖𝑎𝑛 ( 𝑥 𝑖 𝑁-1 , … , 𝑥 𝑗 𝑁-1 ) (Eq. 1)
where x is the number of cells where species i to j that are present in the focal cell occur outside the focal cell, and N is the total number of cells. We subtracted median values from 1 to rank cells with rarer species high and cells with more common species low.

Therefore, SRar is equal to 1 when at least half of the species are only present in the focal cell (i.e., a high degree of rarity), while it is equal to 0 when at least half of the species present in the focal cell also occur in all other cells (i.e., a low degree of rarity). We decided to use the median instead of the mean to reduce the potential effect of outliers, although both measures were highly correlated (Pearson's r = 0.97). The resulting metric was positively related to species richness (R 2 = 0.27; Appendix S7: Fig. S7.1A). We provide an R code to calculate SRar in Appendix S7.

Beta-diversity

We calculated the turnover component of spatial taxonomic (TBDturn) and phylogenetic (PBDturn) beta-diversity of each grid cell following Simpson's index of dissimilarity [START_REF] Baselga | Partitioning the turnover and nestedness components of beta diversity[END_REF][START_REF] Leprieur | Quantifying phylogenetic beta diversity: Distinguishing between 'true' turnover of lineages and phylogenetic diversity gradients[END_REF]:

𝑇𝐵𝐷 𝑡𝑢𝑟𝑛 = 𝑚𝑖𝑛(𝑏,𝑐) 𝑎+𝑚𝑖𝑛(𝑏,𝑐) (Eq. 2)
where a is the number of species shared by both grid cells, b is the number of species exclusive to the focal grid cell, and c is the number of species exclusive to the compared grid cell. The equation for PBDturn is equivalent to that for TBDturn, although shared/exclusive 'species' are replaced by the length of shared/exclusive 'branches' in the phylogenetic tree [START_REF] Leprieur | Quantifying phylogenetic beta diversity: Distinguishing between 'true' turnover of lineages and phylogenetic diversity gradients[END_REF]. The TBDturn and PBDturn metrics reflect true substitution of species or lineages between cells, respectively.

Both the TBDturn and PBDturn metrics have values equal to 0 when two compared grid cells are identical in terms of species and branches, respectively, and values equal to 1 when they are entirely distinct. For each grid cell, turnover was calculated as the mean pairwise dissimilarity between a given cell and all other cells [START_REF] Whittaker | Evolution and measurement of species diversity[END_REF]. This metric is equivalent to the 'local contribution to beta-diversity' proposed by [START_REF] Legendre | Beta diversity as the variance of community data: dissimilarity coefficients and partitioning[END_REF]. We preferred this approach over a moving window approach, which calculates betadiversity among a group of adjacent cells (e.g., [START_REF] Peixoto | Geographical patterns of phylogenetic beta-diversity components in terrestrial mammals[END_REF][START_REF] Qian | Geographic patterns and environmental correlates of taxonomic and phylogenetic beta diversity for large-scale angiosperm assemblages in China[END_REF], because we were interested in comparing areas based on their absolute levels of turnover rather than identifying areas where the local rate of change in turnover was maximum. This corresponds to the difference between the second and the third level of abstraction in betadiversity analyses according to [START_REF] Tuomisto | Analyzing or explaining beta diversity? Understanding the targets of different methods of analysis[END_REF].

We also calculated the deviations of PBDturn given TBDturn (i.e., PBDdev-turn) to identify cells where the phylogenetic turnover was higher or lower than expected given taxonomic turnover [START_REF] Graham | Phylogenetic beta diversity: linking ecological and evolutionary processes across space in time[END_REF][START_REF] Peixoto | Geographical patterns of phylogenetic beta-diversity components in terrestrial mammals[END_REF]:

𝑃𝐵𝐷 𝑑𝑒𝑣-𝑡𝑢𝑟𝑛 = 𝑇𝐵𝐷 𝑡𝑢𝑟𝑛 -(𝑃𝐵𝐷 𝑡𝑢𝑟𝑛 /𝑇𝐵𝐷 𝑡𝑢𝑟𝑛 ) (Eq. 3)
High positive values indicate that phylogenetic turnover is low in relation to taxonomic turnover (i.e., limited lineage exchange), while high negative values indicate that phylogenetic turnover is high in relation to taxonomic turnover (i.e., high lineage exchange).

We calculated all beta-diversity measures using the R package 'betapart' [START_REF] Baselga | betapart: an R package for the study of beta diversity[END_REF].

Statistical analyses

We used structural equation modeling (SEM) to evaluate the effects of environmental and historical variables on TBDturn, PBDturn, and PBDdev-turn. The SEM approach allows causal relationships among variables of interest to be statistically defined and evaluated in the form of mutually interconnected equations [START_REF] Grace | Structural equation modeling and natural systems[END_REF]. In our study system, SEM was deemed as appropriate because of the causal links between beta-diversity and environmental and historical drivers. Specifically, we assumed that these links are mediated by the degree of species uniqueness or rarity (measured here as 'SRar'). Also, we acknowledge that greater environmental heterogeneity, coupled with potential differences in the size of the sampled plots, can influence spatial turnover by affecting regional species richness (SR) [START_REF] Nekola | The distance decay of similarity in biogeography and ecology[END_REF][START_REF] Whittaker | Communities and ecosystems[END_REF]. Therefore, we introduced SR as a predictor of beta-diversity to control for the effect of this variable in our models. Because SR and SRar are nonindependent components of diversity (Appendix S7: Fig. S7.1A), we included a non-causal correlation between SR and SRar in our SEM models.

Our data were spatially structured. We therefore implemented SEM using simultaneous autoregressive models (SAR; [START_REF] Kissling | Spatial autocorrelation and the selection of simultaneous autoregressive models[END_REF], which account for spatial autocorrelation by adding a spatial weight matrix that specifies the neighborhood of each cell and the relative weight of each neighbor. Spatial weight matrices were defined by successively fitting a SAR model and testing several distances between neighbors, ranging from 100 to 500 km away from a given grid cell at intervals of 50 km. We ranked the resulting SAR models by their Akaike information criterion (AIC) values and selected the model with the lowest AIC value as our top SAR model. For all measures, the SAR model with the lowest AIC included neighbors located within 150 km from a given cell. Across all beta-diversity measures, the top SAR model also minimized residual spatial autocorrelation (Appendix S8).

We built a SAR model predicting each of the three beta-diversity metrics (i.e., TBDturn, PBDturn, and PBDdev-turn), based on SR, SRar as well as on the selected environmental (i.e., climate, soil pH, geographical distance from the margin of beech distribution) and historical (i.e., distance from glacial refugia) variables. We modelled betadiversity metrics as a function of the selected environmental and historical variables, rather than as a function of mean geographical and environmental distances because we were interested in exploring how absolute levels of turnover in grid cells varied along biogeographical gradients. We log-transformed the 'distance from the distribution margin', and the 'distance from glacial refugia' to adjust the data distribution closer to normality. We standardized all variables before running the models [START_REF] Zuur | Analysing Ecological Data[END_REF]. To account for the nonlinearity of our predictions, we included quadratic terms, modelled as (x -mean(x)) 2 [START_REF] Maureaud | Biodiversity-ecosystem functioning relationships in fish communities: biomass is related to evenness and the environment, not to species richness[END_REF] for aridity, soil pH and SR. Relationships were included in the final SEM if they were statistically significant and improved the explained variance of each SAR model. The variables included in the models showed low levels of multicollinearity when measured using variance inflation factors (VIF < 2.5; [START_REF] Legendre | Numerical Ecology[END_REF].

To reduce complexity and avoid model saturation of our SEM, we evaluated the goodness-of-fit with Fisher's test after removing non-significant pathways between predictors and beta-diversity metrics. We report Nagelkerke pseudo-R-squared (R 2 ) of final SAR models as a measure of the coefficient of determination. We ran SEM with the R package 'piecewiseSEM' [START_REF] Lefcheck | piecewiseSEM: Piecewise structural equation modeling in R for ecology, evolution, and systematics[END_REF] and SAR models with the R package 'spdep' [START_REF] Bivand | Applied spatial data analysis with R, 2nd edition[END_REF]). We provide an example of the code used for these analyses in Appendix S9. We established significance at P < 0.01 and performed all the analyses in R v. 3.5.3 (R Core Team, 2019).

Results

Effects of environmental and historical factors on taxonomic and phylogenetic turnover

We found the highest values of TBDturn around mountain areas of northern Spain, southern

France, southern Italy, Greece, and more locally in the Netherlands, Belgium, and western Germany (Fig. 1A). The lowest values of TBDturn were found in the northern parts of central Pyrenees, the Jura mountains between France and Switzerland, and in Slovenia. PBDturn showed more homogeneous patterns than TBDturn across Europe, although it highlighted similar regions (Fig. 1B). Taxonomic (TBDturn) and phylogenetic (PBDturn) turnover were highly positively correlated (Pearson's r = 0.78; P < 0.001).

Our final SAR models explained 90% and 63% of the variation in TBDturn and PBDturn, respectively (Fig. 2A andB). The Fisher statistic showed an adequate goodness-offit for both TBDturn (Fisher's C = 12.074; df = 10; P = 0.280) and PBDturn (Fisher's C = 9.203; df = 12; P = 0.685). The TBDturn and PBDturn metrics were best predicted by the species rarity index (SRar) and species richness (SR), which showed a positive linear and a negative quadratic relationship with the two beta-diversity metrics, respectively (Fig. 2A and B; Fig. 3A andB; see also Appendix S10: Fig. S10.2). We also found a negative direct effect of the distance from the distribution margin of beech and the distance from the glacial refugia of beech on TBDturn (Fig. 2B; Fig. 3B; Appendix S11). Soil pH only showed a positive direct effect on PBDturn, and aridity had no significant direct effect on either TBDturn or PBDturn.

When SRar and SR were removed from our models, the explained variation decreased to 46% and 25% for TBDturn and PBDturn, respectively.

The SRar index was negatively related to the distance from the glacial refugia of beech and the distance from the distribution margin of beech but positively related to soil pH (Fig. 2; Appendix S12). The SRar index was lower at intermediate levels of aridity (Fig. 2;

Appendix S12). Species richness was negatively related to aridity and the distance from the glacial refugia of beech but positively related to soil pH (Fig. 2; Appendix S12).

Overall, total effects indicated that both TBDturn and PBDturn were lowest at intermediate levels of aridity but highest towards hotter and drier areas (Fig. 3A andB).

Furthermore, both TBDturn and PBDturn tended to decrease with an increase in soil pH and at a greater distance from the geographical distribution margin of beech. Only TBDturn showed a tendency to decrease at a greater distance from the glacial refugia of beech.

Effects of environmental and historical factors on PBDdev-turn

We observed a strong latitudinal pattern in PBDdev-turn across the geographic distribution of beech (Fig. 1C), as phylogenetic turnover in relation to taxonomic turnover increased towards higher and lower latitudes (see also Appendix S10: Fig. S10.1).

Our final SAR model explained 66% of the variation in PBDdev-turn (Fig. 2C). The Fisher statistic showed an adequate goodness-of-fit (Fisher's C = 7.447; df = 8; P = 0.489).

The PBDdev-turn metric was best predicted by SRar (positively) and SR (negatively) (Fig. 2C; Fig. 3C; see also Appendix S10: Fig. S10.2). We also found a negative direct effect of the distance from the distribution margin of beech, the distance from the glacial refugia of beech, and soil pH on PBDdev-turn (Fig. 2C; Fig. 3C; Appendix S11). All variables had an indirect effect on PBDdev-turn through SRar and SR, although the total effect of aridity was notably lower than for TBDturn and PBDturn. The total effect of the distance from the glacial refugia of beech on PBDdev-turn was notably higher than for TBDturn and showed the opposite pattern than for PBDturn.

Discussion

Our study provides evidence that environmental and historical processes simultaneously drive large-scale patterns of taxonomic and phylogenetic turnover in European beech forests, and that their effects are mediated mainly by the degree of species uniqueness or rarity in regional plant assemblages. Following our initial hypotheses, we demonstrate that (H1) taxonomic and phylogenetic spatial turnover increases in the cells with more distinct environmental conditions or closer to the margin of the geographical distribution of beech, and that (H2) phylogenetic turnover is low in relation to taxonomic turnover closer to the putative glacial refugia of beech, partially because these areas contain a higher fraction of closely related species.

Environmental drivers of taxonomic and phylogenetic turnover

Taxonomic and phylogenetic turnover increased towards both ends of the aridity gradient.

Forests in these conditions are likely to host a relatively high proportion of angiosperm species that can frequently occur in adjacent forest types but only rarely in beech forests. This is in agreement with [START_REF] Hrivnák | Species richness pattern along altitudinal gradient in central European beech forests[END_REF], who showed that the herb-layer of Carpathian beech forests is characterized by two species groups. First, shade-tolerant beech forest specialists are associated with this habitat especially at intermediate altitudes with mild climatic conditions. Second, species typical of neighboring vegetation types are more frequent at both ends of the altitudinal gradient, where the climatic conditions are more extreme. Our study shows that a similar pattern exists at the continental scale and across the full range of European beech forests.

Overall, the highest taxonomic and phylogenetic dissimilarities occurred in hot and dry areas. This high species and lineage turnover might be explained by that fact that such areas generally contain a smaller fraction of the pool of species commonly found in beech forests than those in colder and wetter conditions. The absence of otherwise common species and the presence of uncommon ones, perhaps as a result of strong human influence, might also help explain the high spatial turnover in the relatively species-poor areas in the Netherlands and western Germany (Appendix S13). Our study also indicates that both taxonomic and phylogenetic turnovers increase at intermediate levels of species richness (Appendix S10: Fig. S10.2), suggesting that extremely species-poor beech forests usually host widespread species and lineages.

Beech forests on base-rich soils hosted, on average, more unique sets of species and had a relatively higher taxonomic and phylogenetic spatial turnover. However, the direct positive effect of soil pH on turnover was compensated by an indirect negative effect through species richness. In particular, beech forests on more acidic soils had a lower number of species (see also [START_REF] Ewald | The calcareous riddle: Why are there so many calciphilous species in the Central European flora?[END_REF][START_REF] Jiménez-Alfaro | History and environment shape species pools and community diversity in European beech forests[END_REF][START_REF] Pärtel | Local plant diversity patterns and evolutionary history at the regional scale[END_REF], which in turn promoted a more intense species and lineage turnover than on soils with a higher pH. The strength of these relationships, however, could have been limited by the scale mismatch between the fine-scale effect of soil factors on plants' distribution and the coarse-scale data we used (see also [START_REF] Bruelheide | Global trait-environment relationships of plant communities[END_REF], as well as by the fact that we did not consider variability in soil pH within cells of high geodiversity. Yet, the contrasting effects of soil pH on beta-diversity highlight the need to consider multiple interconnected causal relationships when examining processes structuring natural communities.

As we expected, spatial taxonomic and phylogenetic turnover increased towards the margin of the geographical distribution of beech, reflecting a change in plant composition with proximity to the edges of the habitat range [START_REF] Ries | Ecological responses to habitat edges: Mechanisms, models, and variability explained[END_REF]. Distribution margins frequently represent physical transitions between habitat types and can be characterized by more unusual environmental conditions where edge specialists thrive [START_REF] Laurance | Habitat fragmentation, variable edge effects, and the landscapedivergence hypothesis[END_REF][START_REF] Ries | Ecological responses to habitat edges: Mechanisms, models, and variability explained[END_REF]. However, distribution margins in different locations (e.g., polar, meridional, or continental) can differ fundamentally with respect to historical and other biogeographical factors [START_REF] Bolte | The north-eastern distribution range of European beecha review[END_REF]Weigel et al., 2009). Thus, further research is needed to elucidate how beta-diversity varies among geographical distribution margins in different ecoregions, and whether there is a saturation effect of beta-diversity over distance from these margins.

Both taxonomic and phylogenetic turnovers were structured by similar processes, although with different relative importance. In particular, the explanatory power of our models was higher for taxonomic than for phylogenetic turnover even when we excluded species richness and rarity, indicating that the influence of current climatic conditions on beta-diversity weakens as we go deeper in evolutionary history. A better understanding of the influence of climatic oscillations over geological timescales and of other historical drivers on current patterns of phylogenetic beta-diversity would require complementary biodiversity measures that capture variation in species composition not only at a shallow but also at a deep phylogenetic level (Swenson, 2011a).

Even if based on the largest collection of vegetation-plot data available, our work does not come without uncertainties. First, there is a still a certain amount of unexplained variation in taxonomic and phylogenetic turnover, especially if we remove the effect of species richness and rarity. We suspect this may be due to a combination of local stochastic processes, biotic interactions, distinct dispersal abilities of species and/or unmeasured environmental and spatial variables that ultimately determine the geographical distribution of beech forest species [START_REF] Hubbell | The unified neutral theory of biodiversity and biogeography[END_REF]. Second, although our study covers almost the entire distribution range of European beech forests, further refinement of our conclusions would require additional data in underrepresented areas, such as Germany or western regions of the Balkans (compare Fig. 3 with Appendix S5: Fig. S5.1) and the inclusion of fine-scale ecological factors. Finally, more precise calculations of phylogenetic turnover are needed to challenge our interpretations as molecular data sources provide better and more complete phylogenies.

Historical drivers of taxonomic and phylogenetic turnover

Supporting our initial hypothesis, we found that phylogenetic turnover was less important than taxonomic turnover (i.e., high PBDdev-turn) in the proximity of the glacial refugia of beech and that this effect was largely mediated by species rarity. This finding is consistent with the idea that glacial refugia in mountainous areas host a large fraction of narrow-ranged species [START_REF] Magri | Patterns of post-glacial spread and the extent of glacial refugia of European beech (Fagus sylvatica)[END_REF][START_REF] Willner | Phytogeographical evidence for postglacial dispersal limitation of European beech forest species[END_REF] and that the high replacement of these species is not followed by the replacement of complete phylogenetic lineages.

The observed patterns could likely be explained by the high speciation rate (i.e., vicariance) in mountain refugial areas [START_REF] Graham | Phylogenetic beta diversity: linking ecological and evolutionary processes across space in time[END_REF][START_REF] Turner | Explaining the global biodiversity gradient: energy, area, history and natural selection[END_REF]. For example, repeated glacial isolations have led to recent speciation in some genera like Cardamine or Pulmonaria [START_REF] Cesca | A new species of Cardamine subgen. Dentaria (Cruciferae), apoendemic in Calabria (Southern Italy)[END_REF][START_REF] Kirchner | Molekulare Phylogenie und Biogeographie der Gattung Pulmonaria L. (Boraginaceae)[END_REF], most likely in isolated beech or other deciduous forest environments, resulting in differentiation of the original phylogenetic group into new varieties of species. Furthermore, in Europe, glaciation promoted relatively high extinction rates at northern latitudes [START_REF] Eiserhardt | Climate-driven extinctions shape the phylogenetic structure of temperate tree floras[END_REF][START_REF] Svenning | Could the tree diversity pattern in Europe be generated by postglacial dispersal limitation[END_REF] and allowed the accumulation of more species from ancient lineages in the base-rich mountainous refugial areas of the Iberian, Italian, and Balkan peninsulas [START_REF] Magri | Patterns of post-glacial spread and the extent of glacial refugia of European beech (Fagus sylvatica)[END_REF][START_REF] Willner | Phytogeographical evidence for postglacial dispersal limitation of European beech forest species[END_REF]. Dispersal limitations may have prevented some beech forest specialist species to recolonize areas affected by climate change during interglacial periods [START_REF] Svenning | Could the tree diversity pattern in Europe be generated by postglacial dispersal limitation[END_REF][START_REF] Willner | Phytogeographical evidence for postglacial dispersal limitation of European beech forest species[END_REF]. Altogether, these factors could explain why the putative glacial refugia of beech present a relatively lower variation in phylogenetic turnover in comparison with taxonomic turnover. Nevertheless, we admit that the locations of the glacial refugia of beech used in this study are tentative [START_REF] Magri | Patterns of post-glacial spread and the extent of glacial refugia of European beech (Fagus sylvatica)[END_REF][START_REF] Magri | A new scenario for the Quaternary history of European beech populations: palaeobotanical evidence and genetic consequences[END_REF]. We may be missing other potential refugia, such as cryptic micro-refugia in central Europe [START_REF] Birks | Alpines, trees, and refugia in Europe[END_REF], which may have played a key role in driving current patterns of plant distribution.

Further identification of areas with relatively high rates of speciation (potentially by comparing the coupling between phylogenetic and taxonomic turnover across regions) can help refine the localization of glacial refugia of temperate forests in Europe (see also de [START_REF] De Lafontaine | Beyond skepticism: uncovering cryptic refugia using multiple lines of evidence[END_REF].

Phylogenetic turnover was also less important than taxonomic turnover closer to the glacial refugia of beech after controlling for the effect of species richness and the degree of species rarity in the cells (see also Appendix S11). This suggests that other mechanisms beyond those determining the geographical distribution of beech forest species underly the observed patterns. A plausible explanation is that, once species richness and rarity are held constant, plant assemblages located closer to glacial refugia are phylogenetically more similar to the 'average' assemblage across the range of beech forests. This interpretation is supported by the fact that phylogenetic turnover, unlike taxonomic turnover, decreased with the proximity to glacial refugia. This finding reinforces the idea that glacial refugial areas of European beech forests harbor a large fraction of the plant lineages that can currently be found across beech forests in Europe [START_REF] Petit | Glacial refugia: hotspots but not melting pots of genetic diversity[END_REF].

Conclusions

Previous research on European beech forests has often focused on particular regions or examined the effects of major environmental and historical drivers on alpha-and gammadiversity. Using a large-scale database of vegetation-plot records, we were able to test the effect of environmental and historical evolutionary processes on beta-diversity across the entire distribution range of European beech. We showed that forests located at both ends of the aridity gradient (particularly in warmer and drier areas) and forests located closer to the margins of the geographical distribution of beech hosted more taxonomically and phylogenetically distinct plant assemblages, possibly due to the presence of species associated with other forest types. Moreover, areas located closer to the putative glacial refugia of beech had a comparatively lower phylogenetic turnover in relation to taxonomic turnover. This pattern suggests that these areas still host groups of closely related species that evolved or survived in glacial refugia. We conclude that, despite the vastly coincident geographic patterns in taxonomic and phylogenetic turnover of European beech forests, these two metrics of turnover are distinctly affected by biogeographical processes. Taken together, these findings reinforce the importance of accounting for complementary ecological and evolutionary processes when examining the drivers of spatial beta-diversity at the continental scale.

Biodiversity conservation planning of particular habitats should consider protecting areas at the opposite ends of environmental gradients, as well as those closest to their geographical distribution margin. These areas host unique sets of species that promote greater spatial taxonomic and phylogenetic turnover. Particular attention should be paid to the beech forests located in the warmest and driest areas. These forests have the highest beta-diversity, but their position makes them particularly sensitive to climate change. Further research is thus required to assess the susceptibility of these beech forest communities to climate change and the variation in their distribution, structure, and composition over time. To allow for the comparison of total effects among variables, we placed both positive and negative effects on the same scale and colored them as black and red, respectively. 'Aridity' refers to the first PCA component obtained from the set of climatic variables and depicts an environmental gradient ranging from low temperatures and high precipitation to high temperatures and low precipitation (Appendix S4). The variables 'distance from distribution margin' and the 'distance from glacial refugia' of beech were log-transformed. S4.1). R-squared is also provided. **P < 0.001.

  (see Fig. 9 therein and Appendix S5: Fig. S5.1) and calculated the minimum distance between each cell and these locations to get a measure of distance from refugial areas. Correlations among variables are shown in Appendix S6.

Fig. 1 :

 1 Fig. 1: Spatial patterns of plant beta-diversity in European beech forest vegetation. The maps

Fig. 2 :

 2 Fig. 2: Final structural equation models (SEM) explaining TBDturn (A), PBDturn (B), and

Fig. 3 :

 3 Fig. 3: Total effects of predictor variables on TBDturn (A), PBDturn (B), and PBDdev-turn (C).

Fig

  Fig. S2.1: Histogram of the proportion (%) of non-angiosperm species in the grid cells

Fig. S2. 2 :

 2 Fig. S2.2: Linear relationship between the proportion (%) of non-angiosperm species (log-

Fig. S2. 4 :

 4 Fig. S2.4: Final structural equation models (SEM) explaining TBDturn (A), PBDturn (B), and

Fig. S2. 5 :

 5 Fig. S2.5: Total effects of predictor variables on TBDturn (A), PBDturn (B), and PBDdev-turn (C)

Fig. S3. 1 :

 1 Fig. S3.1: Spatial patterns of beta-diversity of beech forest vegetation in Europe including 10

Fig. S3. 2 :

 2 Fig. S3.2: Total effects of predictor variables on TBDturn (A & D), PBDturn (B & E), and

Fig. S4. 1 :

 1 Fig. S4.1: Map of scores for the first PCA component across the studied grid cells. Areas

Fig. S4. 2 :

 2 Fig. S4.2: Relationship between the three PCA components (PC1, PC2, and PC3) obtained

Fig. S6. 2 :

 2 Fig. S6.2: Linear relationship between the environmental distance of each cell to its centroid

Fig

  Fig. S7.1: Linear relationship between the 'species rarity index' (SRar) and (A) species

Fig. S7. 2 :

 2 Fig. S7.2: Spatial distribution of (A) species richness and (B) the 'species rarity index' (SRar)

Fig. S10. 2 :

 2 Fig. S10.2: Relationship between species richness (SR) (A-C) and the 'species rarity index'

Acknowledgements

We thank Emiliano Agrillo, Henry Brisse, Laura Casella, Panayotis Dimopoulos, Xavier Font, Jean-Claude Gégout, Ute Jandt, Zygmunt Kącki, Remigiusz Pielech, Urban Šilc, Ioannis Tsiripidis, Milan Valachovič, Thomas Wohlgemuth, the Institut national de l'information géographique et forestière (https://inventaire-forestier.ign.fr), and all other database custodians and data contributors for vegetation-plot records; Stephan Hennekens for managing the EVA database; Ilona Knollová for preparing the EVA data for this project, and Irena Axmanová for providing a code to clean up the data.

Funding information

The study was supported by the Czech Science Foundation (project no.19-28491X). I.B. and J.A.C. were funded by the Basque Government (IT936-16).

Data availability statement

The data are available upon request from the European Vegetation Archive (EVA project no. 82).

Author contributions

JPC, BJ-A, FS, WW, and MC conceived the study; IB, JB, JAC, MC, AI, FJ, BJ-A, JL, ZS, KV, and WW collected data; JPC conducted the statistical analyses and wrote the paper; all authors discussed the results and commented on the manuscript.

SUPPORTING INFORMATION

Table of contents:

Appendix S1: Overview of the datasets included in this study.

Table of contents:

Appendix S1: Overview of the datasets included in this study. Appendix S1: Overview of the datasets included in this study. Table S1.1: Overview of the datasets included in this study. For each dataset, we provide the total number and the proportion of plots over the total after resampling procedure (N = 4,930). See the list of databases at http://euroveg.org/eva-database-participating-databases for details. In this study, we obtained, for each plot, data on species occurrences and selected only angiosperms, thus excluding gymnosperms, pteridophytes, lycophytes, and non-vascular plants (~ 4%) to avoid inflated results due to a few taxa connected to deep phylogenetic nodes (e.g., [START_REF] Massante | Environment and evolutionary history depict phylogenetic alpha and beta diversity in the Atlantic coastal white-sand woodlands[END_REF][START_REF] Qian | Geographic patterns and environmental correlates of taxonomic and phylogenetic beta diversity for large-scale angiosperm assemblages in China[END_REF]. Although approximately 87% of the grid cells in our study had less than 15% of non-angiosperm vascular plant species in their species pool (Fig. S2.1), around 10% of the variation in the turnover component of phylogenetic beta-diversity (PBDturn) was explained by the variation in the proportion of nonangiosperm species (log-transformed) in the grid cells (Fig. S2.2).

Database

Our results showed that PBDturn calculated using all vascular plants was less affected by species richness than PBDturn calculated using only angiosperm species (see Fig. S2.5 vs.

Fig. 3). Moreover, soil pH had the opposite effect on PBDturn (i.e., a positive effect) which is largely attributed to the decrease in its indirect effect through species richness.

References: In this section, we present maps and results from SEM models considering thresholds of 10 and 15 plots per grid cell, respectively. Increasing the number of plots per grid cell reduced the total number of grid cells included in our study, but we found no substantial differences from the results produced with 5 per grid cell (see Fig. S3.1 vs. Fig. 1; Fig. S3.2 vs. Fig. 3).

Massante

Appendix S4: Results from Principal Component Analyses (PCA) of environmental variables. Appendix S7: Calculation of the 'species rarity index' (SRar) and relationship with species richness and the 'corrected weighted endemism' (CWE).

The 'species rarity index' (SRar) was positively related to species richness (R 2 = 0.27) and the log-transformed version of the 'corrected weighted endemism' (CWE) (R 2 = 0.43)

proposed by [START_REF] Crisp | Endemism in the Australian flora: Endemism in the Australian flora[END_REF]. The CWE index, as SRar, is also a measure of endemism that is weakly related to species richness. However, CWE emphasizes rarer species over more common species, which is an undesirable property in our study. Therefore, we developed SRar that gives the same weight to all species. We present the correlations between SRar, SR, and the log-transformed version of CWE to clearly show the differences between these metrics.

R code function to calculate the 'species rarity index' (SRar).

data = matrix of species occurrences with sites as rows and species as columns. and the index of species rarity (SRar) derived from the final SEM (Fig. 2). Total effects are defined as the sum of the direct and indirect paths between variables (excluding the paths denoted by double-headed arrows in Fig. 2). To allow for the comparison of total effects among variables, we placed both positive and negative effects on the same scale and colored them as black and red, respectively. 'Aridity' refers to the first PCA component obtained from the set of climatic variables and depicts an environmental gradient ranging from low temperatures and high precipitation to high temperatures and low precipitation (Appendix S4). The variables 'distance from distribution margin' and the 'distance from glacial refugia' of beech were log-transformed.