Karen De Pauw
email: karen.depauw@ugent.be

Camille Meeussen

Sanne Govaert

Pieter Sanczuk

Thomas Vanneste

Markus Bernhardt-Römermann

Kurt Bollmann

Jörg Brunet

Kim Calders

Sara A O Cousins

De Pauw

Edysan Verbeeck

Markus Bernhardt- Römermann

Martin Diekmann

Per- Ola Hedwall

Giovanni Iacopetti

Jonathan Lenoir

Sigrid Lindmo

Anna Orczewska

Quentin Ponette

Jan Plue

Federico Selvi

Fabien Spicher

Hans Verbeeck

Pieter Vermeir

Florian Zellweger

Kris Verheyen

Pieter Vangansbeke

Pieter De Frenne

Taxonomic, phylogenetic and functional diversity of understorey plants respond differently to environmental conditions in European forest edges

Keywords: Taxonomic, phylogenetic and functional diversity of understorey Biodiversity, forest edge, forest understorey, functional diversity, microclimate, phylogenetic diversity, species richness, Determinants of plant community diversity and structure Biodiversity, forest edge, forest understorey, functional diversity, microclimate, phylogenetic diversity, species richness

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Forest biodiversity worldwide is affected by climate change, land-use change and habitat loss [START_REF] Foley | Global consequences of land use[END_REF][START_REF] Lenoir | A significant upward shift in plant species optimum elevation during the 20th century[END_REF][START_REF] Pereira | Global Biodiversity Change: The Bad, the Good, and the Unknown[END_REF], Vellend et al. 2013, Zellweger et al. 2020). Future climate change is predicted to cause further biodiversity losses [START_REF] Thomas | Extinction risk from climate change[END_REF][START_REF] Malcolm | Global warming and extinctions of endemic species from biodiversity hotspots[END_REF][START_REF] Thuiller | Consequences of climate change on the tree of life in Europe[END_REF][START_REF] Trisos | The projected timing of abrupt ecological disruption from climate change[END_REF] as biodiversity redistribution is hampered by habitat fragmentation in terrestrial systems [START_REF] Lenoir | Species better track climate warming in the oceans than on land[END_REF]. In forests, the largest part of plant species richness, up to more than 80%, is located in the understorey [START_REF] Gilliam | The ecological significance of the herbaceous layer in temperate forest ecosystems[END_REF]. Furthermore, understorey communities play an important role in forest ecosystem dynamics by mediating nutrient cycling, tree regeneration and other crucial ecosystem functions [START_REF] Gilliam | The ecological significance of the herbaceous layer in temperate forest ecosystems[END_REF], Landuyt et al. 2019). Resurvey studies showed no general decline in species richness of understorey communities over the past decades due to balanced local colonisations and extinctions [START_REF] Keith | Taxonomic homogenization of woodland plant communities over 70 years[END_REF], Vellend et al. 2013). However, if a limited set of taxa replaces many different species across the forest biome, homogenization can lead to losses of total biodiversity on the scale of the forest biome, even though locally no decline in species richness is registered [START_REF] Staude | Replacements of small-by large-ranged species scale up to diversity loss in Europe's temperate forest biome[END_REF]. During the past century, herbaceous understorey plant communities homogenized by an increasing presence of nutrient-demanding and shade-tolerant species [START_REF] Keith | Taxonomic homogenization of woodland plant communities over 70 years[END_REF][START_REF] Naaf | Habitat specialists and generalists drive homogenization and differentiation of temperate forest plant communities at the regional scale[END_REF][START_REF] Prach | Landscape-scale vegetation homogenization in Central European submontane forests over the past 50years[END_REF][START_REF] Van Den Berge | Contrasting vegetation change (1974-2015) in hedgerows and forests in an intensively used agricultural landscape[END_REF][START_REF] Staude | Replacements of small-by large-ranged species scale up to diversity loss in Europe's temperate forest biome[END_REF]. Simultaneously, climate change causes an increasing dominance of warm-adapted species in understorey communities, a process referred to as thermophilization [START_REF] Bertrand | Changes in plant community composition lag behind climate warming in lowland forests[END_REF], De Frenne et al. 2013, Zellweger et al. 2020). Such changes in community composition might affect the functional or phylogenetic diversity of understorey communities, as well as the role of the understorey in forest ecosystem dynamics [START_REF] Wardle | Terrestrial Ecosystem Responses to Species Gains and Losses[END_REF]. Therefore, it is important to expand our understanding from species richness to functional and phylogenetic diversity and to investigate which environmental factors are driving different diversity patterns in understorey plant communities. Human pressure on forests leads to forest fragmentation, and consequently to increasing forest edge to interior ratios with important consequences for forest biodiversity [START_REF] Haddad | Habitat fragmentation and its lasting impact on Earth's ecosystems[END_REF] and its redistribution as climate warms [START_REF] Lenoir | Species better track climate warming in the oceans than on land[END_REF]. Currently, 70% of forested area is located closer than one km to a forest edge and 20% is even closer than 100 m [START_REF] Haddad | Habitat fragmentation and its lasting impact on Earth's ecosystems[END_REF]. Therefore, it is important to understand the ecological processes occurring in forest edges, in addition to those of forest interiors. Many environmental factors change drastically from the forest's edge to its interior [START_REF] Matlack | Microenvironment variation within and among forest edge sites in the eastern United-States[END_REF][START_REF] Gehlhausen | Vegetation and microclimatic edge effects in two mixed-mesophytic forest fragments[END_REF]. Forest edges receive more atmospheric acidifying and eutrophying deposition [START_REF] Devlaeminck | Variation in throughfall deposition across a deciduous beech (Fagus sylvatica L.) forest edge in Flanders[END_REF][START_REF] Wuyts | Comparison of forest edge effects on throughfall deposition in different forest types[END_REF]) and have higher nitrogen and carbon stocks, compared to the interior [START_REF] Remy | Strong gradients in nitrogen and carbon stocks at temperate forest edges[END_REF][START_REF] Meeussen | Drivers of carbon stocks in forest edges across Europe[END_REF]. Therefore, the increasing proportion of forest edges might accelerate the current homogenization of understorey plant communities, which is characterised by an increasing presence of nutrient-demanding species [START_REF] Verheyen | Driving factors behind the eutrophication signal in understorey plant communities of deciduous temperate forests[END_REF][START_REF] Van Den Berge | Contrasting vegetation change (1974-2015) in hedgerows and forests in an intensively used agricultural landscape[END_REF][START_REF] Staude | Replacements of small-by large-ranged species scale up to diversity loss in Europe's temperate forest biome[END_REF]. Furthermore, forest edges are characterised by increased wind speeds and incoming solar radiation, resulting in more variable microclimates and drier soil conditions than forest interiors [START_REF] Matlack | Microenvironment variation within and among forest edge sites in the eastern United-States[END_REF][START_REF] Chen | Microclimate in forest ecosystem and landscape ecology -Variations in local climate can be used to monitor and compare the effects of different management regimes[END_REF][START_REF] Gehlhausen | Vegetation and microclimatic edge effects in two mixed-mesophytic forest fragments[END_REF]. Such conditions typically harbour communities with a high proportion of generalist species, as forest generalists avoid the shady, humid and strongly buffered microclimate of forest interiors [START_REF] Normann | How forest edge-center transitions in the herb layer interact with beech dominance versus tree diversity[END_REF], Govaert et al. 2020). Such edge conditions and a dominance of forest generalists might locally reduce the functional or phylogenetic diversity of these forest edge communities. However, it is not yet clear which environmental factors affect functional and phylogenetic diversity of understories in forest edges. The macroclimate and surrounding landscape matrix have an additional effect on understorey diversity next to edaphic and stand conditions. Macroclimate is a well-known driver of biodiversity [START_REF] Francis | A globally consistent richness-climate relationship for angiosperms[END_REF]Currie 2003, Kreft and[START_REF] Kreft | Global patterns and determinants of vascular plant diversity[END_REF]. Mean annual temperature and potential evapotranspiration, for example, are important predictors of plant species richness [START_REF] Francis | A globally consistent richness-climate relationship for angiosperms[END_REF][START_REF] Kreft | Global patterns and determinants of vascular plant diversity[END_REF][START_REF] Qian | A latitudinal gradient in large-scale beta diversity for vascular plants in North America[END_REF]. In Europe, many forests are situated in fragmented landscapes and many forest specialists have low colonization capacities [START_REF] Verheyen | Response of forest plant species to land-use change: a life-history trait-based approach[END_REF][START_REF] Hermy | Legacies of the past in the present-day forest biodiversity: a review of past land-use effects on forest plant species composition and diversity[END_REF]; the amount of habitat and fragmentation can thus have a considerable effect on understorey diversity [START_REF] Valdes | The contribution of patch-scale conditions is greater than that of macroclimate in explaining local plant diversity in fragmented forests across Europe[END_REF], Govaert et al. 2020). The habitat amount hypothesis states that not necessarily patch isolation or size, but the amount of habitat present in the 'local landscape' affects species density [START_REF] Fahrig | Rethinking patch size and isolation effects: the habitat amount hypothesis[END_REF][START_REF] Watling | Support for the habitat amount hypothesis from a global synthesis of species density studies[END_REF]. Lower species richness of forest generalists and especially specialists were found in forest patches surrounded by less forested area (within a radius of 100-500 m) [START_REF] Valdes | The contribution of patch-scale conditions is greater than that of macroclimate in explaining local plant diversity in fragmented forests across Europe[END_REF][START_REF] Takkis | Drivers of species richness and community integrity of small forest patches in an agricultural landscape[END_REF], Govaert et al. 2020). Most studies on understorey biodiversity rely on taxonomic species richness due to its simplicity and convenience. However, during the last two decades, the focus has changed from number of species towards their ecological diversity, i.e. the degree to which species differ in terms of their function, niche or evolutionary history [START_REF] Cadotte | The ecology of differences: assessing community assembly with trait and evolutionary distances[END_REF]). The ecological diversity can be assessed with a functional approach, based on functional traits, and with a phylogenetic approach, based on species' genealogies. Both approaches can add complementary information for conservation (Carvalho et al. 2017, Cadotte and[START_REF] Cadotte | Difficult decisions: Strategies for conservation prioritization when taxonomic, phylogenetic and functional diversity are not spatially congruent[END_REF]. The functional diversity metrics can provide important information regarding ecosystem functioning [START_REF] Cadotte | Beyond species: functional diversity and the maintenance of ecological processes and services[END_REF][START_REF] Flynn | Functional and phylogenetic diversity as predictors of biodiversity-ecosystem-function relationships[END_REF], whereas phylogenetic diversity relates to genetic variability, which is believed to improve the communities' adaptability to future change [START_REF] Cavender-Bares | The merging of community ecology and phylogenetic biology[END_REF]. Functional and phylogenetic metrics can provide insight in the community assembly of forest understorey communities [START_REF] Gerhold | Functional and phylogenetic community assembly linked to changes in species diversity in a long-term resource manipulation experiment[END_REF][START_REF] Thorn | Changes in the dominant assembly mechanism drive species loss caused by declining resources[END_REF][START_REF] Vanneste | Functional trait variation of forest understorey plant communities across Europe[END_REF]. At the finest spatial resolution of forest plant communities (usually 400 m 2 ), community assembly is often attributed to the limiting similarity and competitive exclusion mechanisms [START_REF] Webb | Phylogenies and community ecology[END_REF][START_REF] Cavender-Bares | The merging of community ecology and phylogenetic biology[END_REF]). These mechanisms suggest that closely related (e.g. sister species) or functionally similar species compete more intensely than phylogenetically or functionally distant species. Consequently, the chance of co-existence is higher for distantrelated or functionally divergent species [START_REF] Webb | Phylogenies and community ecology[END_REF]. On a larger scale involving landscape or regional extents and coarser spatial resolutions (between 400 m 2 and 1 km 2 ), the environmental filtering mechanism limits the diversity of communities, filtering for species adapted to the specific environment through similar ecological strategies and/or phylogenetic histories [START_REF] Cavender-Bares | The merging of community ecology and phylogenetic biology[END_REF][START_REF] Laliberte | Environmental filtering explains variation in plant diversity along resource gradients[END_REF]. When species retain their niche and related ecological traits over time, with niche defined as the set of biotic and abiotic conditions where a species can persist [START_REF] Holt | Bringing the Hutchinsonian niche into the 21st century: Ecological and evolutionary perspectives[END_REF]), this can be described as niche conservatism [START_REF] Wiens | Niche conservatism as an emerging principle in ecology and conservation biology[END_REF]. Phylogenetic niche conservatism then expands this concept to related species [START_REF] Wiens | Niche conservatism as an emerging principle in ecology and conservation biology[END_REF]. The degree of phylogenetic conservatism of functional traits and the shape of the phylogenetic tree thus determine the correspondence of phylogenetic and functional diversity approaches [START_REF] Mazel | Conserving Phylogenetic Diversity Can Be a Poor Strategy for Conserving Functional Diversity[END_REF]. They might result in contrasting plant diversity patterns and therefore, contribute different important information to provide a more general view on the communities' biodiversity [START_REF] Cadotte | The ecology of differences: assessing community assembly with trait and evolutionary distances[END_REF][START_REF] Thorn | Rare species, functional groups, and evolutionary lineages drive successional trajectories in disturbed forests[END_REF]. While numerous studies have unravelled how taxonomic and compositional diversity are affected in deciduous forests by environmental drivers [START_REF] Van Calster | Diverging effects of overstorey conversion scenarios on the understorey vegetation in a former coppice-with-standards forest[END_REF][START_REF] Price | Small-scale patterns of species richness and floristic composition in relation to microsite variation in herb-rich woodlands[END_REF], Depauw et al. 2019, Macek et al. 2019), only few studies have simultaneously assessed taxonomic, functional and phylogenetic diversity of the understorey (e.g. [START_REF] Wasof | Dominance of individual plant species is more important than diversity in explaining plant biomass in the forest understorey[END_REF] and [START_REF] Closset-Kopp | Do drivers of forestry vehicles also drive herb layer changes (1970-2015) in a temperate forest with contrasting habitat and management conditions?[END_REF]) and never before at the continental scale. Here we assessed different metrics of taxonomic, phylogenetic and functional diversity of understorey herb communities in 225 plots spread along edge-to-interior gradients in deciduous forests across Europe. We capitalized on large-scale macroclimatic latitudinal and altitudinal gradients in combination with fine-grain management and edge-to-core gradients in microclimate. We specifically assessed the following hypotheses for the understorey biodiversity in European forest edges:

H1: We expected different responses to environmental predictors, in terms of magnitude and direction, between the taxonomic, phylogenetic and functional diversity metrics that we derived in our study. H2: We expected taxonomic diversity to increase with higher light availability and less buffered microclimates, in contrast to functional and phylogenetic diversity, which we expected to decrease due to the increased presence and dominance of generalist species: i.e. local functional and phylogenetic homogenization. H3: We expected a higher importance of edaphic and stand conditions than landscape conditions for phylogenetic and functional diversity, as edaphic and stand conditions strongly influence community assembly processes such as environmental filtering and species' competitive ability.

Methods

Study area and experimental set-up

We selected European broadleaved forests with a dominance of oak species (mainly Quercus robur, Q. petraea and Q. cerris) locally complemented with Fagus sylvatica, Betula pubescens, Populus tremula, Ulmus glabra, Alnus incana and Carpinus betulus, because these are important for conservation as biodiversity 'hotspots'.

All forests had a minimal area of 4 ha and were ancient forests; they have been continuously forested and were not converted to other land use since at least the oldest available maps, which is typically at least 150-300 years. We specifically did not include post-agricultural forests, to rule out any possible effect of past land-use history. [START_REF] Karger | Data Descriptor: Climatologies at high resolution for the earth's land surface areas[END_REF]). The latitudinal range of the broadleaved temperate forest biome, as given by [START_REF] Olson | Terrestrial Ecoregions of the World: A New Map of Life on Earth: A new global map of terrestrial ecoregions provides an innovative tool for conserving biodiversity[END_REF], was covered completely. In three of these regions (Italy, Belgium and Norway), an additional elevational gradient was established; forests were selected at low, intermediate and high elevations (ranging from 21 to 908 m above sea level, corresponding to a smaller macroclimate gradient of 1.5 to 4 °C MAT). In the other six regions, only lowland forests were selected with elevations between 7.5 and 451 m above sea level (Fig.

SI.D1).

In each of these regions, and at each elevation, forest stands with three different management types were selected. The first type were 'dense forests', which were not thinned for at least 10 to 30 years. Additional criteria for this forest type were a well-developed shrub layer and a complex vertical structure. These forests generally had low canopy openness (5.8 ± 0.6%, mean of three densiometer measurements [START_REF] Baudry | Estimating light climate in forest with the convex densiometer: operator effect, geometry and relation to diffuse light[END_REF]) and high basal area (mean and standard error of 28.8 ± 1.5 m²/ha). The second forest type, referred to as 'intermediate forests', were forests that were regularly thinned with the most recent thinning ideally 5 to 10 years ago. For this forest type, we looked for sparser shrub layers and a less complex vertical structure (canopy openness of 6.5 ± 0.6% and mean basal area of 31.4 ± 1.9 m²/ha). The third forest type were 'open forests', which were regularly thinned and most recently within 4 years before sampling. Additionally, the shrub layer and subdominant tree layer of these forests were ideally sparse or lacking and the vertical structure comprised only the dominant tree layer. These forests were generally characterised by high canopy openness values (14.8 ± 2.1%) and low basal area values (21.6 ± 1.3 m²/ha).

In each forest stand, a transect was established perpendicular to the forest edge. Forest edge was defined as the outer edge of the forest stand that borders a matrix of open land. All forest edges were south-facing to ensure comparability since edge orientation can affect forest microclimate and herbaceous vegetation considerably [START_REF] Matlack | Microenvironment variation within and among forest edge sites in the eastern United-States[END_REF][START_REF] Honnay | Permeability of ancient forest edges for weedy plant species invasion[END_REF][START_REF] Orczewska | Floristic analysis of the two woodland-meadow ecotones differing in orientation of the forest edge[END_REF]. Five plots of 3 by 3 m were installed, with their plot centre on an exponentially increasing distance from the forest edge, respectively 1.5, 4.5, 12.5, 35.5 and 99.5 m. All plots were at least 100 m away from any forest edge other than the studied forest edge, to avoid interference with effects from other forest edges. Thus, in total 45 forest edge-interior transects were sampled ((6 lowland regions + 9 from 3 regions with elevation gradient) x 3 forest types), totalling 225 plots (45 forest edges x 5 plots per edge-interior transect). More detailed descriptions of the selection criteria and structural characteristics of the different management types can be found in [START_REF] Govaert | Edge influence on understorey plant communities depends on forest management[END_REF] and [START_REF] Meeussen | Structural variation of forest edges across Europe[END_REF].

Trait and phylogeny data Vegetation surveys were performed during peak of the vegetation season (May -July 2018) according to the local phenology. In each 3 x 3 m plot, all vascular plant species were identified and their percentage ground cover was estimated (n = 353 species in total). The herb layer comprised all vascular plants smaller than 1m, including woody, non-woody plants and lianas. Vegetation surveys are also described in detail by Govaert et al. (2020). Seedlings, shrub species and lianas in the herb layer were excluded from this analysis (n = 62 species), since they do not remain in the herb layer throughout their lifecycle and their trait values from most online databases do not represent the juveniles encountered in the herb layer (Table SI.D1).

Three key functional traits were chosen based on the leaf-height-seed plant ecology strategy scheme: seed mass, specific leaf area (SLA) and plant height. SLA informs on the plant's acquisition and conservation of resources as it represents a trade-off between photosynthetic rate and leaf lifespan [START_REF] Wright | The worldwide leaf economics spectrum[END_REF]). Plant height strongly determines the plant's competitiveness for light and the dispersal of its seeds, whereas seed mass reflects the trade-off between seed output and seedling survival [START_REF] Westoby | A leaf-height-seed (LHS) plant ecology strategy scheme[END_REF][START_REF] Westoby | Land-plant ecology on the basis of functional traits[END_REF][START_REF] Diaz | The global spectrum of plant form and function[END_REF]. We assessed the variation of traits in our study species with a principal component analysis (PCA). A parallel analysis for determining significant principal components (method developed by Franklin et al. (1995)) showed that two significant principal components are necessary to represent the variation in the trait data (Fig. SI.D2a). The first two axes comprise 40.9 and 34.4% of the trait variation and the biplot shows the separate factor loading of functional traits on these axes (Fig. SI.D2b). The plot shows a broad spread of species, representing different plant strategies in three major trait domains informing on species' resource use, competition and reproduction [START_REF] Westoby | A leaf-height-seed (LHS) plant ecology strategy scheme[END_REF][START_REF] Westoby | Plant ecological strategies: Some leading dimensions of variation between species[END_REF][START_REF] Pierce | How well do seed production traits correlate with leaf traits, whole-plant traits and plant ecological strategies?[END_REF], Diaz et al. 2016) (Fig. SI.D2b). Species-specific trait values were derived from several databases including the LEDA trait database (Kleyer et al. 2008), BiolFlor (Kuhn et al. 2004) and the Kew Seed Information Database (KEW 2017) (Table SI.D4). For 60 species, no or insufficient trait data was present, these species were excluded from the analysis (Table SI.D4). However, trait values were available for all species that occurred in more than 5% of the plots with a mean cover value of more than 2% (Table SI.D2).

We chose to extract a phylogenetic tree from the dated molecular phylogeny for land plants constructed by Zanne et al. (2014), because this tree included plants from the entire study region, including Mediterranean species. The extraction was done with the 'brranching' package in R [START_REF] Chamberlain | brranching: Fetch 'Phylogenies' from Many Sources[END_REF]. Two species were omitted from further analysis due to their absence in the phylogenetic tree of Zanne et al. (2014) (Phegopteris connectilis and Polypodium interjectum). Both species occurred in less than 3% of the plots. The resulting final vegetation matrix of the herb layer contained 229 species (Fig. 1), which represented on average 93.7% of the total herbaceous cover in the plots (Table SI.D3). The final tree counted 229 tips, one for each species in the final vegetation matrix, and had 221 internal nodes (Fig. 1). Figure 1. Phylogenetic tree of species included in this study. This tree was standardised taxonomically with The Plant List (2013) and visualized with ggtree and gheatmap in R [START_REF] Yu | GGTREE: an R package for visualization and annotation of phylogenetic trees with their covariates and other associated data[END_REF]. Functional traits are represented with a colour scale around the tree, plant height (m) on the inner circle, SLA (mm 2 mg -1 ) on the middle circle and seed mass (mg) on the outer circle.

Diversity metrics

Traditionally, the diversity concept includes two components: species richness and evenness. Species richness gives the number of species, and species evenness the equitability of their relative abundances. Generally, species diversity is defined as a metric including both species richness and species evenness [START_REF] Smith | A consumer's guide to evenness indices[END_REF]. However, based on phylogeny or function, diversity can also be regarded as a measure of variation in the community, irrespective of species richness. Therefore, we always calculated three metrics, i.e., richness, evenness and diversity or variability metrics, and this based on taxonomy, phylogeny or function to come to nine response variables (Table 1).

Regarding taxonomy, we calculated species richness (Tax.rich) (number of species per plot) and the Evar evenness index (Tax.even), proposed by [START_REF] Smith | A consumer's guide to evenness indices[END_REF]. This evenness index ranges from zero to one, is independent from species richness and was calculated with the 'codyn' R package [START_REF] Hallett | CODYN: AnR package of community dynamics metrics[END_REF].

As diversity metric we calculated the Shannon diversity index (Tax.div) [START_REF] Shannon | A Mathematical Theory of Communication[END_REF]. Phylogenetic diversity metrics were calculated, using the 'pez' R package [START_REF] Pearse | pez: phylogenetics for the environmental sciences[END_REF]. First, the phylogenetic species variability (Phy.div) of the community indicates the variation in evolutionary history of the species. The metric ranges between zero and one and is independent of the community's species richness.

A value of one characterizes a community in which none of the species has a lineage in common, and a value close to zero indicates a community of species, which share large parts of their lineages in the phylogenetic tree [START_REF] Helmus | Phylogenetic measures of biodiversity[END_REF]). Secondly, phylogenetic species richness (Phy.rich) was calculated by multiplying Phy.div with species richness values [START_REF] Helmus | Phylogenetic measures of biodiversity[END_REF]). Thirdly, phylogenetic species evenness (Phy.even) was calculated by adapting Phy.div to take the abundances of species into account, as such it combines evenness in abundance and phylogeny [START_REF] Helmus | Phylogenetic measures of biodiversity[END_REF]). Finally, functional metrics were calculated using the 'FD' package in R [START_REF] Laliberté | A distance-based framework for measuring functional diversity from multiple traits[END_REF]. Traits were standardized to mean zero and unit variance and a species-species Euclidian distance matrix was computed. The 'Cailliez' correction method was used to correct for negative eigenvalues [START_REF] Cailliez | The analytical solution of the additive constant problem[END_REF]) and a Principal Coordinates Analysis (PCA) was performed with the resulting species-species distance matrix. The axes obtained from the PCA were used to compute functional richness (Fun.rich) and functional evenness (Fun.even) indices [START_REF] Villéger | New multidimensional functional diversity indices for a multifaceted framework in functional ecology[END_REF]. Functional richness reports on the trait-space volume, whereas functional evenness assesses simultaneously the evenness of species distribution in trait-space and evenness of their abundances. Finally, the Rao's quadratic entropy (Fun.div) was calculated based on relative species abundances and pairwise functional differences between species (Botta-Dukat 2005, Laliberté and Legendre 2010). This distance-based metric is frequently used to quantify functional diversity and is independent of species richness (Botta-Dukat 2005, Laliberté and Legendre 2010). Table 1. Overview of the diversity metrics used to assess the diversity of understorey plants in 225 plots in the forest edges of European deciduous forests. A richness, evenness and diversity or variability metric (.rich, .even, .div respectively) were calculated based on taxonomy, phylogeny or function (Tax., Phy., Fun. respectively). In the formulas, p i , p j = relative abundance of species i, j, S = total number of species in community and d ij = difference between species i, j. The value of the Pearson correlation is given as a number and indicated by the colour in the Pearson correlation matrix. The mean pairwise functional difference between species in the community, weighted by their abundance.

𝐹𝑢𝑛.𝑑𝑖𝑣 = 𝑆 -1 ∑ 𝑖 𝑆 ∑ 𝑗

Environmental predictor variables Edaphic conditions

In each plot of 3 × 3 m, the forest floor or organic soil horizon (i.e. litter, humus and fragmentation layer) was sampled in a 20 × 20 cm subplot from its surface to the mineral soil layer underneath, after removal of the herb layer. These samples were dried to constant weight at 65 °C for 48 h and then, the organic soil layer mass was determined (mass OS). This variable was used as indicator for litter quality, thickness and nutrient availability since low-degradable litter tends to accumulate on the forest floor and results in slower nutrient turnover and lower nutrient availability [START_REF] Scott | Foliage litter quality and annual net N mineralization: Comparison across North American forest sites[END_REF]. Additionally, dense litter layers may pose a physical barrier for germination of forest species or reduce germination through phytotoxic components [START_REF] Facelli | Plant litter-it's dynamics and effects on plant community structure[END_REF]. Texture analysis (% sand, silt and clay) was performed by sieving and sedimentation of mineral soil samples of 10-20 cm depth, whereas soil pH was determined for mineral soil samples of 0-10 cm depth. For these analyses, five subsamples were taken per plot and pooled (detailed description in Suppl. Information A).

Stand conditions

In each plot, the microclimate temperature was recorded hourly at 1 m height using a temperature data logger (Lascar EL-USB-1, range of -30 to +80 °C, resolution of 0.5 °C) covered by a radiation shield (Fig.

SI.D3).

For each of the nine regions and for each of the three elevation levels, the temperature was also measured in an identical set-up in an open field close to the forest stands ('reference' sensor). The temperature measurements of these 'reference' sensors were used to calculate temperature offset values (offset = subcanopy temperature -free air temperature = plot sensor -'reference' sensor). Positive and negative offset values represent warmer and cooler forest microclimates, respectively, compared to the macroclimate temperature. Typically, the forest microclimate is buffered from temperature extremes and this buffering is largest during the summer months, when understorey plants are most likely to experience extreme heat and drought stress [START_REF] Zellweger | Seasonal drivers of understorey temperature buffering in temperate deciduous forests across Europe[END_REF]. Also cold temperatures and frost are limiting factors regarding plant survival and distributions [START_REF] Sakai | Frost survival of plants: Responses and Adaptation to Freezing Stress[END_REF][START_REF] Woodward | The impact of low temperatures in controlling the geographical distribution of plants[END_REF][START_REF] Svenning | Postglacial dispersal limitation of widespread forest plant species in nemoral Europe[END_REF][START_REF] Bucher | Temporal and spatial trade-offs between resistance and performance traits in herbaceous plant species[END_REF].

For these reasons, we focused on the effect of cooling of maximum temperatures in summer and warming of minimum temperatures in winter. During winter (from October 2018 to March 2019), the offset was calculated for the mean daily 5 th percentile temperature ('winter offset') and during summer (from April to September 2019) for the mean daily 95 th percentile temperatures ('summer offset'). The three forest management types are expected to impact microclimate and light availability at the forest floor, since they differ in density and complexity. We quantified forest structural differences using Plant Area Index (PAI), which is half of the surface area of all aboveground vegetation matter (including stems, branches and leaves) per unit surface area. PAI was calculated as the integral of vertically resolved plant area per volume density (m² m -3 ) profiles derived from single-scan position terrestrial laser scanning (TLS) using a RIEGL VZ-400 (RIEGL Laser Measurement Systems GmbH, Horn, Austria) in the centre of each plot. PAI thus gives an indication of the denseness and complexity of the forest structure and is negatively related to light availability at the forest floor. We chose to apply TLS, since the technique is highly reproducible and more direct compared to conventional forest surveys [START_REF] Calders | Monitoring spring phenology with high temporal resolution terrestrial LiDAR measurements[END_REF][START_REF] Liang | Terrestrial laser scanning in forest inventories[END_REF][START_REF] Calders | Variability and bias in active and passive ground-based measurements of effective plant, wood and leaf area index[END_REF]). The TLS method was described in detail by [START_REF] Meeussen | Structural variation of forest edges across Europe[END_REF].

Landscape conditions

The macroclimate was taken into account as mean annual temperature and annual precipitation, which were retrieved for the coordinates of each plot from the CHELSA database for 1979 -2013 (resolution of ~ 1 km 2 ) [START_REF] Karger | Data Descriptor: Climatologies at high resolution for the earth's land surface areas[END_REF]). Plots from one edge-to-interior transect might be located within the same or neighbouring 1 km 2 grid cells and thus have very similar macroclimatic conditions, whereas microclimatic conditions will vary depending on the forest edge distance and structure. To incorporate the amount of habitat in the 'local landscape', the percentage area with a tree cover >20% was calculated within a radius of 500 m based on satellite-based global tree cover data with a spatial resolution of 30 m [START_REF] Hansen | High-Resolution Global Maps of 21st-Century Forest Cover Change[END_REF].

Mean and standard deviations of the nine predictor variables (mass OS, sand fraction, pH, plant area index, winter offset, summer offset, mean annual temperature, annual precipitation, forest cover) are given for the nine study regions in Table SI.D5.

Data analyses

We used generalized linear mixed-effects models (GLMM) to infer responses of the nine diversity metrics (Table 1) to environmental drivers and performed all analyses in R (R Development Core Team 2020). Due to the hierarchical nature of the data, GLMMs with transect ID (45 levels corresponding to 45 edge-to-interior transects) as random effect (random intercept) nested within region (nine levels corresponding to 9 regions) were used (225 plots nested in 45 transect nested in 9 regions). For models with species richness as response variable, a Poisson error distribution was used with a log link function since these are count data, and, as a consequence, these models can be nonlinear. For all other models a Gaussian error distribution was applied, resulting in strictly linear models. Correlations between predictor variables were assessed with Pearson correlation coefficients before modelling (Fig. SI.D4). Multicollinearity of the predictor variables in the models was assessed using variance inflation factors (VIFs) with the vif function from the package 'cars' [START_REF] Fox | An R Companion to Applied Regression[END_REF]. For all models, VIFs were smaller than 3.1 and thus no strong multi-collinearity issues were detected among the set of predictor variables we used [START_REF] Neter | Applied linear statistical models. Regression, analysis of variance, and experimental design[END_REF][START_REF] Zuur | Mixed effects models and extensions in ecology with R[END_REF]). The models were fitted with the 'lme4' package [START_REF] Bates | Fitting Linear Mixed-Effects Models Using lme4[END_REF].

To explore the effects of latitude, elevation, forest structure type and the distance to the edge (the design variables of the experimental set up) on the diversity metrics, GLMM were performed with these design variables as fixed effects, including all two-way interactions. The distance to the edge was log-transformed to meet model assumptions (log e ). Furthermore, we included continuous environmental variables in the models. To represent the edaphic conditions, we included sand fraction, pH of the mineral soil and mass of the organic soil layer (mass OS).

Mass OS was log-transformed due to its skewed distribution (log e ). For forest structure, we used the PAI (Plant Area Index), summer offset and winter offset. For landscape conditions, we included two macroclimate variables, mean annual temperature (MAT) and annual precipitation, and percentage forest cover. No interactions were taken into account to avoid too much predictor terms and complexity. Equation 1 represents the model structure of the global model.

y ~ % 𝑠𝑎𝑛𝑑 + 𝑝𝐻 + log (𝑚𝑎𝑠𝑠 𝑂𝑆) + 𝑠𝑢𝑚𝑚𝑒𝑟 𝑜𝑓𝑓𝑠𝑒𝑡 + 𝑤𝑖𝑛𝑡𝑒𝑟 𝑜𝑓𝑓𝑠𝑒𝑡 + 𝑃𝐴𝐼 + 𝑀𝐴𝑇 + 𝑎𝑛𝑛𝑢𝑎𝑙 𝑝𝑟𝑒𝑐𝑖𝑝𝑖𝑡𝑎𝑡𝑖𝑜𝑛 + % 𝑓𝑜𝑟𝑒𝑠𝑡 𝑐𝑜𝑣𝑒𝑟 + (1|𝑟𝑒𝑔𝑖𝑜𝑛/𝑡𝑟𝑎𝑛𝑠𝑒𝑐𝑡)
Eq. 1 Starting from the global model with all predictor variables, model selection was performed based on the lowest corrected Akaike Information Criterion (AICc), testing all possible combinations of predictor variables with the dredge function from the package 'MuMin' [START_REF] Barton | MuMIn: Multi-Model Inference[END_REF]. During model selection, maximum likelihood was used to fit models, whereas afterwards, the restricted maximum likelihood approach was used to obtain model estimates of the best model. P-values were obtained with the 'lmerTest' package [START_REF] Kuznetsova | lmerTest Package: Tests in Linear Mixed Effects Models[END_REF] and corrected for multiple testing with false discovery rates based on Pike (2011). All continuous predictor variables were scaled to unit variance and mean zero. As second step in the analysis, we performed a variation partitioning to assess the importance of the predictor variables of edaphic, stand and landscape conditions, following the procedure of [START_REF] Legendre | Numerical ecology[END_REF]. This second step was only performed for Tax.rich, Phy.div and Fun.div, as we chose to focus further analysis on this set of three independent metrics to obtain comparable and independent results regarding taxonomic, phylogenetic and functional diversity (Table 1) [START_REF] Botta-Dukat | Rao's quadratic entropy as a measure of functional diversity based on multiple traits[END_REF][START_REF] Helmus | Phylogenetic measures of biodiversity[END_REF][START_REF] Schleuter | A user's guide to functional diversity indices[END_REF]. First, the global model was produced, including all nine predictors (Eq. 1), and then models were produced with only one or two out of three groups (edaphic -stand -landscape conditions). Each group had three predictor variables to balance the variation partitioning. For all models we obtained the proportion of variance explained by the fixed factors of the model (marginal R 2 , R 2 m ), according to Nakagawa and Schielzeth (2013). Finally, we calculated the amount of variation explained by each group and combination of groups by subtracting R 2 m from the R 2 m of the global model. The variation explained by each group and intersection is reported as percentage of variation explained by the global model.

Results

The results that we highlight here are strongly focused on taxonomic richness (Tax.rich), phylogenetic diversity (Phy.div) and functional diversity (Fun.div) because they are most widely used and independent from changes in species number [START_REF] Botta-Dukat | Rao's quadratic entropy as a measure of functional diversity based on multiple traits[END_REF][START_REF] Helmus | Phylogenetic measures of biodiversity[END_REF][START_REF] Schleuter | A user's guide to functional diversity indices[END_REF]).

Diversity patterns with latitude, forest type and distance to forest edge We found that latitude, forest type and the distance to the forest edge, as well as the interaction between forest type and the edge distance, significantly explained variation in several biodiversity metrics (Table SI.B1). Tax.rich decreased towards the forest interior, especially strong close to the edge (estimate ± standard error: -0.096 ± 0.022, p < 0.001). Contrastingly, Phy.div decreased towards the forest edge and this gradient was the steepest for the intermediate forest type, whereas dense forests had a higher Phy.div overall and the smallest edge-to-interior gradient (interaction estimate ± se: -0.045 ± 0.017, p = 0.029). While Tax.rich and Phy.div showed no latitudinal gradient, Fun.div decreased towards the north (estimate ± se: -0.005 ± 0.001, p = 0.005). Simultaneously, there was a significant difference between the forest management types, with the intermediate type showing the highest Fun.div (open compared to intermediate type estimate ± se: -0.009 ± 0.003, p = 0.03) (Fig. 2). Figure 2. Taxonomic richness, phylogenetic diversity and functional diversity as a function of the latitude (a-c), distance to the forest edge (d-f) and the interaction with forest management type (c, e). The lines show model predictions for significant effects based on the generalized mixed effect models (Table SI.B1) and shading corresponds to 95% confidence intervals. Yellow, blue and red colours in d and e indicate significant (interactive) effects of forest management type (see legend). Jittering on the X axis was added for clarity, as well as transparency of points, darker areas thus indicate several overlapping points. Model fit is shown in the panels with significant predictor as marginal R 2 (R 2 m ) and conditional R 2 (R 2 c ), following Nakagawa and Schielzeth (2013). We highlight results for Tax.rich combined with Phy.div and Fun.div because of their wide use and independence from species number [START_REF] Botta-Dukat | Rao's quadratic entropy as a measure of functional diversity based on multiple traits[END_REF][START_REF] Helmus | Phylogenetic measures of biodiversity[END_REF][START_REF] Schleuter | A user's guide to functional diversity indices[END_REF]).

Diversity patterns with landscape, stand and edaphic environmental conditions

In general, we found four out of nine environmental predictors having a significant impact on multiple diversity metrics: soil pH, Plant Area Index (PAI), summer offset and mean annual temperature (MAT). From these, only summer offset explained variation for taxonomic, phylogenetic and functional diversity metrics. Annual precipitation and winter offset were not retained for any of the biodiversity metrics after model selection. The percentage forest cover in the surrounding landscape, mass of the organic soil layer and sand fraction showed no significant effects after correcting for multiple testing (Table 2). Furthermore, none of the environmental predictors significantly explained variation in functional richness and functional evenness, neither for taxonomic and phylogenetic evenness after p-value correction (Table 2). Tax.rich decreased with decreasing light availability (increasing PAI) and stronger microclimatic buffering of summer maximum temperature (more negative summer offset), whereas Fun.div increased when microclimatic buffering was stronger (Fig. 3 b and h, Table 2). Additionally, Tax.rich increased for higher pH values, whereas Phy.div decreased with increasing pH values (Fig. 3 c and f, Table 2). Furthermore, Fun.div increased with increasing MAT, while Tax.rich and Phy.div showed no significant response to MAT (Fig. 3 a, d and g, Table 2). Tax.rich, Phy.div and Fun.div thus exhibited contrasting responses to stand and edaphic conditions (Fig. 3). Figure 3. The relationship of taxonomic richness, phylogenetic and functional diversity with landscape, stand and edaphic conditions as predictors (Table 2). Taxonomic richness (Tax.rich), phylogenetic diversity (Phy.div) and functional diversity (Fun.div) as a function of mean annual temperature (a,d,g), summer offset (b,e,h) and soil pH (c,f,i). The lines show model predictions for significant parameter estimates based on the generalized mixed effect models and shading corresponds to 95% confidence intervals. Jittering was added for clarity, as well as transparency of points, darker areas thus indicate several accumulated points at the same or overlapping location. We highlight results for Tax.rich combined with Phy.div and Fun.div because of their wide use and independence from species number [START_REF] Botta-Dukat | Rao's quadratic entropy as a measure of functional diversity based on multiple traits[END_REF][START_REF] Helmus | Phylogenetic measures of biodiversity[END_REF][START_REF] Schleuter | A user's guide to functional diversity indices[END_REF]).

Variation partitioning: landscape, stand and edaphic conditions For Tax.rich, all three variable groups (landscape, stand and edaphic conditions) explained large parts of the variation (42.8%, 52.4% and 34.1% respectively). For both Phy.div and Fun.div we did find large differences between proportions of variation explained by the different groups. Most of the variation in Phy.div was explained by edaphic conditions (52.6%), followed by stand conditions (24.1%) and only a small amount of variation was explained by the landscape conditions (8.5%). For Fun.div, the variation partitioning resulted in an opposite pattern. Landscape conditions explained most of the variation in Fun.div (60.2%), followed by the stand conditions (23.6%) and only a small amount of variation was explained by edaphic conditions (3.2%) (Fig. 4, Fig. SI.D5). Figure 4. Variation partitioning. Venn-Euler diagrams for variation partitioning of taxonomic richness, phylogenetic and functional diversity. These diagrams show the proportion of explained variation (marginal R 2 ) by three variable groups (landscape, stand and edaphic conditions) and the shared proportion of group combinations (intersection of circles) compared to the explained variation of the global model including all nine predictors. The size and numbers within the circles correspond to the proportion of explained variation. The proportion shared by all three groups is indicated with an arrow. Table 2. Summary of the results of the generalized linear mixed-effect models of the nine diversity metrics with landscape, stand and edaphic conditions as predictors. Models were run for taxonomic, phylogenetic and functional richness (Tax.rich, Phy.rich, Fun.rich), evenness (Tax.even, Phy.even, Fun.even) and diversity metrics (Tax.div, Phy.div, Fun.div). Model fit was assessed based on marginal R 2 (R 2 m ), the proportion of variation explained by the fixed effects, and conditional R 2 (R 2 c ), the proportion of variation explained by both random and fixed effects (Nakagawa and Schielzeth 2013). Parameter estimates are given for the model with lowest corrected AIC value after model selection with the dredge function. P-values corrected for false discovery rates are given following Pike (2011), those estimates with a p-value below 0.05 are given in bold. The significance of the original p-values is given within brackets with '***' for p <0.001, '**' for p <0.01, '*' for p <0.05. Positive and negative significant parameter estimates are respectively red and blue coloured, with darker hues for significant terms after correction.

Predictor

Landscape

Discussion

Responses differ between taxonomic, phylogenetic and functional diversity metrics (H1)

The use of phylogenetic and functional diversity metrics along taxonomic metrics provided different, complementary information on community assembly mechanisms in forest understorey plant communities across Europe. From the Pearson correlation matrix, it was clear that some metrics were strongly correlated, such as taxonomic and phylogenetic richness, but that others were not (Table 1). Both phylogenetic and functional diversity were independent of taxonomic richness, showing that in this study, plots with higher species numbers in the understorey did not necessarily have higher phylogenetic or functional understorey diversity. Furthermore, we found different responses in both magnitude and direction to gradients of the experimental design of the study (Fig. 2, Table SI.B1) and to other environmental predictors (Table 2). Low mean annual temperatures (e.g. at high latitudes) led to lower functional diversity, whereas taxonomic richness and phylogenetic diversity showed no clear pattern driven by the macroclimate (Table 2, Fig. 2, Fig. 3). At high latitudes, environmental filtering probably led to a plant community with restricted variation in functional traits, enabling the plants to survive longer and colder winters. For example, the community weighted mean of plant height tended to be smaller and SLA tended to be higher in plots with lower mean annual temperatures. Additionally, colder, more northern regions needed to be recolonized after the last glaciation period, causing a convergent functional composition regarding functional traits linked to dispersal, such as lower seed mass values [START_REF] Pinto-Ledezma | Patterns of Beta Diversity of Vascular Plants and Their Correspondence With Biome Boundaries Across North America[END_REF]. Also soil pH caused contrasting responses; species richness decreased with a lower pH, possibly due to higher soluble aluminium in acidic soils, which decreases plant root growth [START_REF] Kopittke | Identification of the Primary Lesion of Toxic Aluminum in Plant Roots[END_REF][START_REF] Bojórquez-Quintal | Aluminum, a Friend or Foe of Higher Plants in Acid Soils[END_REF], in combination with the toxic effect of a low pH itself [START_REF] Falkengrengrerup | Experimental evidence for the relative sensitivity of deciduous forest plants to high soil acidity[END_REF]Tyler 1993, Falkengren-Grerup 1995). Phylogenetic diversity, on the other hand, increased when the soil became more acidic, corroborating the findings of [START_REF] Piwczynski | Influence of tree plantations on the phylogenetic structure of understorey plant communities[END_REF], which might be due to environmental filtering for acidophilous species. Most often, environmental filtering leads to phylogenetic clustering, as was found in the understorey of Mediterranean oak forests [START_REF] Selvi | Impact of pine invasion on the taxonomic and phylogenetic diversity of a relict Mediterranean forest ecosystem[END_REF]). This pattern can arise when the ability to cope with a particular environmental filter is shared in more closely related species of the community. In our study, however, these acidophilous communities often contained a mixture of forbs, graminoids and ferns (for example, Stellaria holostea and Pteridium aquilinum co-occurring with Molinia caerulea in plots with low pH). Species from different phylogenetic lineages acquired the ability to grow in acidic soils independently, leading to phylogenetic variability instead. Clearly, the effect of environmental filtering on phylogenetic diversity depends strongly on the degree of phylogenetic conservatism of the traits which are favoured by the environmental filter. These results suggest that in this situation niche conservatism might be the case within species but not phylogenetically, among related species.

There was no clear effect of soil pH on functional diversity. Plants might cope with acidity in different ways, which are not necessarily directly related to their general ecological strategy in terms of resource use, competition and reproduction (represented by SLA, plant height and seed mass). Mechanisms to deal with acidity are often related to root traits, such as the secretion of organic acid anions or associations with mycorrhiza [START_REF] Marschner | Mechanisms of adaptation of plants to acid soils[END_REF][START_REF] Chen | Adaptation to acidic soil is achieved by increased numbers of cis-acting elements regulating ALMT1 expression in Holcus lanatus[END_REF]. Furthermore, root traits provide highly useful information on, for example, velocity of resource turn-over and the association with mycorrhiza [START_REF] Bergmann | The fungal collaboration gradient dominates the root economics space in plants[END_REF].

Whereas some categorical root traits are becoming increasingly available (such as mycorrhiza type), quantitative root traits are not yet available for many understorey herbs. For example, specific root length data was only available for 37% of the species studied here, whereas this is one of the most studied belowground traits (Table SI.D7). In the future, root traits could greatly increase our understanding of biodiversity responses to environmental conditions, especially to edaphic factors such as soil pH.

Effect of light availability and microclimatic buffering (H2) Higher taxonomic richness was found in plots with higher light availability, a pattern clearly driving the increase towards the forest edges through a gradual decrease in plant area index (Table 2, Fig. 2, Fig. SI.D6) [START_REF] Gehlhausen | Vegetation and microclimatic edge effects in two mixed-mesophytic forest fragments[END_REF][START_REF] Honnay | Permeability of ancient forest edges for weedy plant species invasion[END_REF][START_REF] Vallet | The effects of urban or rural landscape context and distance from the edge on native woodland plant communities[END_REF]). The higher understorey species richness at the forest edge in this study is predominantly driven by an increase of forest generalist species [START_REF] Govaert | Edge influence on understorey plant communities depends on forest management[END_REF]. Furthermore, the generalist richness was lower beneath canopy tree species casting more shade [START_REF] Govaert | Edge influence on understorey plant communities depends on forest management[END_REF]. The lower light availability in the forest interior poses an environmental filter for the light demanding generalist edge species [START_REF] Pellissier | Understorey plant species show long-range spatial patterns in forest patches according to distanceto-edge[END_REF]. The increasing availability of a limiting resource, such as light, can remove an environmental filter, increasing the number of species, but this process can lead to different edge-to-interior patterns depending on the context and land-use history (e.g. [START_REF] Hofmeister | Strong influence of long-distance edge effect on herb-layer vegetation in forest fragments in an agricultural landscape[END_REF]).

We also detect higher species richness in the less microclimatically buffered forest edge plots (Fig. 2, Fig.

SI.D7

). Strikingly, the functional diversity was lower in plots that are less buffered (Fig. 3). Similarly, the plots in the open forest management type had lower functional diversity and less buffered summer temperatures (Fig. 2, Fig. SI.D7). It is possible that reduced microclimatic buffering due to less dense tree canopies made the forest understorey more susceptible to spring frost and summer drought [START_REF] Von Arx | Microclimate in forests with varying leaf area index and soil moisture: potential implications for seedling establishment in a changing climate[END_REF], Zellweger et al. 2019), acting as an environmental filter and favouring plants adapted to higher temperatures and lower soil moisture. Additionally, it is known that thermophilization of understorey communities (the process in which warm-adapted species gain abundance over more cold-adapted species over time) occurs more in forest stands with less microclimatic buffering, especially when light availability is high (De Frenne et al. 2013).

Thermophilization can be driven by tall and competitive species [START_REF] De Frenne | Light accelerates plant responses to warming[END_REF], which could suppress smaller, less competitive plants and drive a loss of functional diversity through competitive exclusion. To assess this hypothesis, resurveyed vegetation studies are needed (De Frenne et al. 2013[START_REF] Feeley | Climate-driven changes in the composition of New World plant communities[END_REF], Zellweger et al. 2020). Nevertheless, in an additional analysis we calculated the floristic temperature of the understorey communities and assessed its relationship with microclimate. The floristic temperature increased when minimum winter temperatures were more buffered (less cold in winter), but we did not detect a link with the buffering of summer maximum temperatures (SI.C).

Besides, forest microclimates are more buffered in warmer regions (Fig. SI.D7, De Frenne et al. ( 2019)), and functional diversity increased with warmer macroclimate temperatures; this association could thus have contributed indirectly to the higher functional diversity in strongly buffered forest plots. Conversely, unmanaged forests often have a high basal area, high PAI and a complex forest structure providing strongly buffered microclimates [START_REF] Frey | Spatial models reveal the microclimatic buffering capacity of old-growth forests[END_REF]) and low light levels on the forest floor [START_REF] Hardwick | The relationship between leaf area index and microclimate in tropical forest and oil palm plantation: Forest disturbance drives changes in microclimate[END_REF]. Higher functional diversity is often found in dense, unmanaged forests with low light levels [START_REF] Liu | Understorey plant community and light availability in conifer plantations and natural hardwood forests in Taiwan[END_REF], Closset-Kopp et al. 2019[START_REF] Lelli | Biodiversity response to forest structure and management: Comparing species richness, conservation relevant species and functional diversity as metrics in forest conservation[END_REF]). In such conditions, the proportion of generalist species is lower [START_REF] Govaert | Edge influence on understorey plant communities depends on forest management[END_REF]) and the less competitive forest specialists shape the understorey communities [START_REF] Honnay | Permeability of ancient forest edges for weedy plant species invasion[END_REF]).

Chesson's framework states that coexistence mechanisms can result from stabilizing or equalizing processes. species traits; but see HilleRisLambers et al. ( 2012)) can be ensured, if biotic and abiotic factors force species to experience stronger intraspecific than interspecific competition [START_REF] Chesson | Mechanisms of maintenance of species diversity[END_REF]. Equalizing niche differences, on the other hand, are processes that suppress fitness inequality between species resulting from many ecological and evolutionary factors, which fundamentally contributes to stable multispecies coexistence. Low light availability can function as an equalizing mechanism, decreasing the growth of the most competitive species. Most forest specialist species might be quite similar in having low competitive abilities, e.g. lower plant height than forest generalist species [START_REF] Marinšek | What makes a plant species specialist in mixed broadleaved deciduous forests?[END_REF]. As a result, the interspecific competition can be low and thus small niche differences could already cause larger intraspecific competition and stabilize the coexistence between these species. This is supported by the tendency of evenness metrics to increase with lower light availability. Furthermore, stronger microclimatic buffering is often related to a higher complexity of the forest structure [START_REF] Frey | Spatial models reveal the microclimatic buffering capacity of old-growth forests[END_REF][START_REF] Kovacs | Stand structural drivers of microclimate in mature temperate mixed forests[END_REF], which enhances heterogeneity in forest stratification and in forest-floor conditions. This could also increase the coexistence of species with different abiotic preferences and explain the increased functional diversity. Further research might look deeper into the link between microclimatic buffering, thermophilization, functional diversity in forests and possible consequences for ecosystem functioning, such as nutrient cycling and litter decomposition.

We hypothesised that increased light and reduced microclimatic buffering could lead to local functional or phylogenetic homogenization due to the increased presence of generalist species. Contrary to our expectation, we did not find lower functional diversity in forest edges or with more light (Fig.SI.B1, Table 2). However, weaker buffering of the microclimate did decrease functional diversity (Fig. 3, Table 2). Further research might elucidate the different roles of light and microclimate in shaping the diversity of the forest understorey through controlled experimental set-ups separating both factors.

Regarding the phylogenetic diversity, we detected a significant interaction of the forest management type and distance to the edge (Fig. 2). For the intermediate forest type, the phylogenetic diversity decreased towards the forest edge, whereas it was more stable for the dense and open forest type. Both the increase of generalist species and the decrease of Plant Area Index at the forest edge were steepest and most abrupt for the intermediate forest types [START_REF] Govaert | Edge influence on understorey plant communities depends on forest management[END_REF], Meeussen et al. 2020). Abrupt changes in forest structure could thus be related to changes in phylogenetic diversity in the understorey. Stand characteristics did account for 24% of explained variation of phylogenetic diversity in the variation partitioning, however, no strong direct linear responses to microclimate or plant area index were found after model selection (Fig. 4, Table 2).

Varying importance of landscape, stand and edaphic conditions (H3) Landscape, stand and edaphic conditions were of relatively similar importance for Tax.rich. This could be expected from the well-known influence of edaphic and stand conditions on species richness [START_REF] Van Calster | Diverging effects of overstorey conversion scenarios on the understorey vegetation in a former coppice-with-standards forest[END_REF][START_REF] Vanhellemont | Relating changes in understorey diversity to environmental drivers in an ancient forest in northern Belgium[END_REF], Govaert et al. 2020) in combination with the large spatial gradient of the study, covering important changes in landscape conditions [START_REF] Bernhardt-Römermann | Drivers of temporal changes in temperate forest plant diversity vary across spatial scales[END_REF]. The phylogenetic and functional diversity metrics showed very different results. For Phy.div, the importance of the landscape conditions was negligible compared to the influence of stand and edaphic conditions. However, [START_REF] Li | Species richness and phylogenetic diversity of native and nonnative species respond differently to area and environmental factors[END_REF] reported a high impact of macroclimatic factors for phylogenetic diversity but, in contrast to our study, their study extended to all ecosystem types and covered a larger spatial gradient, with larger variation in macroclimate. For the functional diversity, the macroclimatic temperature gradient explained most of the variation, acting as a strong overarching environmental filter, followed by the forest microclimate. We expected also a strong impact of edaphic conditions, similar to phylogenetic diversity, but this was not confirmed. The recurring difference between functional and phylogenetic diversity responses indicates that the conservation of plant height, SLA and seed mass through the phylogenetic tree of the understorey species is not very strong (Fig. 1), challenging the niche conservatism concept and the notion of functional and phylogenetic diversity as substitutes. Furthermore, it also shows that the evolutionary history indeed comprises much more information about a species than three key functional traits, even though they represent the main plant strategies regarding resource acquisition, competitiveness and reproduction. Stand conditions were important for the three diversity metrics through regulation of light at the forest floor and microclimatic buffering. Nevertheless, the directions of responses were different. Trade-offs between diversity aspects are clearly present in deciduous forests. It is thus essential to look beyond species counts when studying forest understorey plant diversity, certainly when interested in community assembly mechanisms.

Conclusions

Combining taxonomy, phylogeny and functional traits proved to be important when assessing plant diversity of understorey communities in forest edges of temperate deciduous forests. We show that different diversity aspects can be driven by contrasting environmental conditions and in different directions, leading to trade-offs between diversity metrics. It was clear that also for forest understorey species, the reality is more complex than environmental filtering and competitive exclusion leading unequivocally to less and more diverse communities respectively. Diversity is an outcome of many different, interacting and context-dependent processes. Functional and phylogenetic diversity were no mere substitutes for each other, but revealed complementary and important insights. Future studies could acknowledge this complexity more by including intraspecific trait variation, an aspect which we did not consider here, but could help to understand functional biodiversity patterns and community assembly mechanisms [START_REF] Siefert | Incorporating intraspecific variation in tests of trait-based community assembly[END_REF][START_REF] Violle | The return of the variance: intraspecific variability in community ecology[END_REF][START_REF] Des Roches | The ecological importance of intraspecific variation[END_REF]. Furthermore, we suggest to measure and study belowground traits of understorey herbs to enable a clearer understanding of functional diversity responses to edaphic conditions. We detected no significant local functional or phylogenetic homogenization close to forest edges in general. However, we did find a decrease of phylogenetic diversity in forest edges for forests of intermediate and high openness and we did find a decreasing functional diversity in plots with less buffered microclimates. In the context of climate change, with increasing frequency of extreme summer temperatures (IPCC 2018) and canopy disturbances due to drought, heat stress and insect attacks [START_REF] Allen | A global overview of drought and heat-induced tree mortality reveals emerging climate change risks for forests[END_REF][START_REF] Anderegg | Tree mortality from drought, insects, and their interactions in a changing climate[END_REF], it would be highly valuable to further investigate this relationship and possible consequences of functional homogenization in the understorey for forest ecosystem functioning. ), the proportion of variation explained by the fixed effects, and conditional R 2 (R 2 c ), the proportion of variation explained by both random and fixed effects (Nakagawa and Schielzeth 2013). Parameter estimates are given for the model with lowest corrected AIC value after model selection with the dredge function. Estimates with a p-value below 0.05 are given in bold. Positive and negative significant parameter estimates are respectively red and blue coloured, with darker hues for significant terms.

Figuresb

  Figures

  Figure SI.D6 Plant Area Index (PAI) given for the forest type and edge-to-interior gradients used in

  deviation of the nine continuous predictor variables for the different regions. IT = Italy, SW = Switzerland, FR = France, BE = Belgium, PO = Poland, GE = Germany, SS = Southern Sweden, CS = Central Sweden, NO = Norway. MAT = Mean annual temperature (°C), AP = annual precipitation (mm), PAI = Plant Area Index, mass(OS) = mass of organic soil layer (g).

  To increase comparability, all forests had loamy soils with an intermediate moisture content.

We selected forests in nine regions spanning a 2000-km long latitudinal gradient across Europe (from south to north): Italy, Switzerland, France, Belgium, Poland, Germany, southern Sweden, central Sweden and Norway. This latitudinal gradient includes a change in mean annual temperature (MAT) of > 10°C and a variation in annual precipitation from approximately 550 to 1250 mm (data retrieved from CHELSA database for 1979 -2013 (resolution of ~ 1 km 2 )

Table SI .B1. Summary of the results of the linear mixed effect models of the nine diversity metrics with latitude, elevation, forest management type and distance to forest edge as predictors.

 SI Intermediate forests were used as the reference forest management type. Model fit was assessed based on marginal R 2 (R 2 m ), the proportion of variation explained by the fixed effects, and conditional R 2 (R 2 c ), the proportion of variation explained by both random and fixed effects(Nakagawa and Schielzeth 2013). Parameter estimates are given for the model with lowest AIC value after model selection with the dredge function. P-values corrected for false discovery rates are given following(Pike 2011) , those estimates with a p-value below 0.05 are given in bold. The significance of the original p-values is additionally given within brackets with '***' for p <0.001, '**' for p <0.01, '*' for p <0.05. Positive parameter estimates are red coloured, negative parameter estimates are blue coloured, non-significant parameters are not coloured. Lat = latitude, Dist= distance to forest edge.

	Response		Latitude	Elevation	Log(Distance) Open type	Dense type	Dist * Lat	Dist * Open	Dist * Dense R 2 m	R 2 c
	Tax.rich	Mean± SE			-0.096±0.022					0.02	0.77
		p corrected (original)			<0.001 (***)
	Tax.even	Mean ± SE p corrected (original)		-0.014±0.02 0.592						0.004 0.77	0.35
	Tax.div	Mean ± SE			-0.066±0.027					0.01	0.68
		p corrected (original)			0.047 (*)
	Phy.rich	Mean ± SE	1.129±0.291							0.00 0.21	0.76
		p corrected (original)	0.005 (***)
	Phy.even	Mean ± SE p corrected (original)	0.027±0.013 0.073 (*)		0.035±0.010 0.005 (***)			0.023±0.009 0.047 (*)		0.10 0.22	0.22
	Phy.div	Mean ± SE			0.057±0.012	0.024±0.038	0.052±0.038		-0.031±0.018 -0.045±0.017 0.09	0.50
		p corrected (original)			<0.001 (***)	0.537	0.198		0.103	0.029 (**)
	Fun.rich	Mean ± SE			-0.050±0.016 -0.052±0.041 -0.028±0.041		0.061±0.023 0.053±0.023	0.01 0.04	0.51
		p corrected (original)			0.012(**)	0.250	0.514		0.027 (**)	0.046 (*)
	Fun.even	Mean ± SE p corrected (original)								0 0.68	0.26
	Fun.div	Mean ± SE	-0.005±0.001			-0.009±0.003 -0.006±0.003			0.23	0.61
		p corrected (original)	0.005 (***)			0.030 (**)	0.430
										0.21
										0.76

Table SI.C Summary of the results of the generalized linear mixed-effect model of the floristic temperature with landscape, stand and edaphic conditions as predictors.

  Model fit was assessed based on marginal R 2 (R 2 m

Table SI .D2 Herb species excluded from the analysis due to lack of trait data or phylogenetic data. Only 101

 SI two species were not present in the phylogenetic tree ofZanne et al. (2014), these are marked with *. Other 102 species were excluded due to lack of sufficient trait data. The number and percentage of plots they were present 103 in and the mean and maximum cover they took in these plots are given.

		Ranunculus velutinus	7	0.031	0.01	0.01
		Rosa sempervirens	8	0.036	0.01	0.01
		Rubus canescens	1	0.004	0.01	0.01
		Rubus hirtus	12	0.053	0.02	0.06
		Rubia peregrina	9	0.040	0.01	0.04
		Saussurea alpina	3	0.013	0.01	0.03
		Secale cereale	1	0.004	0.01	0.01
		Excluded herb species Sherardia arvensis Smyrnium perfoliatum Symphytum bulbosum	Number plots 4 1 8	% of plots 0.018 0.004 0.036	mean % cover 0.01 0.02 0.01	0.01 max % 0.02 cover 0.01
	104	Alchemilla alpina Thesium linophyllon Allium pendulinum Trifolium ochroleucum Anemone apennina Urospermum dalechampii Angelica archangelica Veronica agrestis Aristolochia rotunda Veronica arvensis Asparagus acutifolius Veronica persica Asplenium adiantum.nigrum Vicia bithynica Asparagus tenuifolius Vicia ochroleuca Bistorta vivipara Viola alba Campanula rapunculus Viola biflora Carex distachya Phegopteris connectilis * Carex olbiensis Polypodium interjectum * Carex ovalis mean of all excluded herb species Cirsium heterophyllum max of all excluded herb species Cistus salviifolius	3 1 7 4 10 4 2 1 2 1 8 2 3 1 1 6 1 28 3 8 5 5 10 3 1 4.71 1 28 5	0.013 0.004 0.031 0.018 0.044 0.018 0.009 0.004 0.009 0.004 0.036 0.009 0.013 0.004 0.004 0.027 0.004 0.124 0.013 0.036 0.022 0.022 0.044 0.013 0.004 0.021 0.004 0.124 0.022	0.01 0.01 0.01 0.01 0.10 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.02 0.01 0.01 0.01 0.01 0.04 0.01 0.12 0.02 0.01 0.03 0.01 0.01 0.12 0.05	0.01 0.02 0.01 0.01 0.02 0.30 0.01 0.01 0.01 0.01 0.01 0.05 0.01 0.01 0.01 0.02 0.07 0.01 0.10 0.01 0.20 0.04 0.01 0.07 0.03 0.03 0.30 0.01 0.10
	105	Colchicum neapolitanum	1	0.004	0.01	0.01
	106	Crepis leontodontoides Crocus etruscus	8 7	0.036 0.031	0.01 0.01	0.02 0.01
		Crocus vernus	8	0.036	0.01	0.01
		Cyclamen hederifolium	3	0.013	0.01	0.01
		Cyclamen repandum	17	0.076	0.01	0.05
		Daphne gnidium	1	0.004	0.01	0.01
		Dorycnium hirsutum	1	0.004	0.01	0.01
		Echinops ritro	1	0.004	0.01	0.01
		Empetrum nigrum	1	0.004	0.01	0.01
		Geranium dissectum	1	0.004	0.01	0.01
		Geranium purpureum	11	0.049	0.01	0.01
		Geranium rotundifolium	1	0.004	0.01	0.01
		Gnaphalium uliginosum	2	0.009	0.01	0.01
		Helianthemum nummularium	1	0.004	0.01	0.01
		Lathyrus venetus	8	0.036	0.04	0.25
		Lilium bulbiferum	3	0.013	0.01	0.01
		Lolium perenne	3	0.013	0.01	0.01
		Oenanthe pimpinelloides	15	0.067	0.02	0.10
		Phleum pratense	4	0.018	0.01	0.03
		Polygala flavescens	4	0.018	0.02	0.03
		Potentilla micrantha	5	0.022	0.01	0.01
		Prunella laciniata	2	0.009	0.01	0.01
		Pulmonaria hirta	1	0.004	0.01	0.01
		Pulicaria odora	3	0.013	0.01	0.01

Table SI .D3 Percentage groundcover of the total herbaceous layer from species included in the 107 analysis per region.

 SI

	Region	Mean % cover
		of plots
	Norway	0.993
	Central Sweden	0.982
	Southern Sweden	1
	Germany	1
	Belgium	0.997
	Poland	1
	France	0.999
	Switzerland	1
	Italy	0.707
	108
	109

Table SI.D4 Source databases for trait data.

  BiolFlor:(Kuhn et al. 2004); D-3 The Dispersal and Diaspore Database[START_REF] Hintze | D-3: The Dispersal and Diaspore Database -Baseline data and statistics on seed dispersal[END_REF]; The Ecological Flora database[START_REF] Fitter | The ecological flora database[END_REF]; LEDA trait database:(Kleyer et al. 2008); O: Own measurements; Rothmaler, Exkursionsflora von Deutschland[START_REF] Jäger | Rothmaler, Exkursionsflora von Deutschland[END_REF]; SID: Seed Information Database Kew: (KEW 2017). Trait values were available for all species that occurred in more than 5% of the plots with a mean cover value of more than 2%. LEDA, BiolFlor, SID, Ecological Flora TableSI.

	Trait	Units	Data source
	Plant height	m	LEDA, Ecological Flora, Rothmaler, O
	Specific leaf area (SLA) mm 2 /mg LEDA, O
	Seed mass	mg	D-3,

D5 Mean and standard

Table SI .D6 Mean and standard deviation of the nine response variables for the different regions

 SI . IT = Italy, SW = Switzerland, FR = France, BE = Belgium, 123 PO = Poland, GE = Germany, SS = Southern Sweden, CS = Central Sweden, NO = Norway. Fine.root.length.per.fine.root.dry.mass.. specific.fine.root.length..SRL.. Root.tissue.density..root.dry.mass. per.root.volume.. Fine.root.nitrogen..N..content. per.fine.root.dry.mass. Fine.root.tissue.density.. fine.root.dry.mass.per.fine.root.volume.. Fine.root.dry.mass.per.fine.root.fresh.mass.. fine.root.dry.matter.content..RDMC.. Plant.biomass.and.allometry.. Root.mass.fraction..fraction.of.root.dry.mass. per.whole.plant.dry.mass..

	122	.Fine.root.diameter.		TRY		26.64	61
		.TRY		26.64	61
		Root_C_concentration		GRooT	26.64	61
		Tax.rich .Root.diameter.	Tax.even	TRY	Tax.div 26.20	60	Phy.rich	Phy.even	Phy.div
		MAT (°C) 10.84 ± 3.44 .TRY AP (mm) forest cover (%) IT 0.79 ± 0.12 1.84 ± 0.37 25.76 SW 3.60 ± 1.76 0.63 ± 0.29 0.82 ± 0.42 IT 13.03 ± 1.78 758.22 ± 104.27 69.33 ± 14.54 FR 5.93 ± 2.91 0.54 ± 0.29 0.93 ± 0.62 Root_C_N_ratio GRooT 21.40	summer offset (°C) 4.04 ± 1.30 59 1.54 ± 0.54 -6.99 ± 2.94 2.78 ± 1.58 49	winter offset 0.32 ± 0.05 0.24 ± 0.14 (°C) 1.96 ± 2.16 0.28 ± 0.16	0.37 ± 0.03 PAI 0.39 ± 0.06 3.79 ± 0.81 0.46 ± 0.14
		SW BE .Plant.biomass.and.allometry.. 10.46 ± 0.08 910.13 ± 4.55 5.25 ± 4.28 0.66 ± 0.26	70.02 ± 10.73 1.03 ± 0.64 TRY 20.09	-8.49 ± 0.95 2.52 ± 1.46 46	0.51 ± 0.12 0.35 ± 0.22	5.09 ± 1.21 0.52 ± 0.22
		FR PO Root.dry.mass.per.plant. 10.71 ± 0.23 10.80 ± 4.07	681.47 ± 44.89 0.65 ± 0.13	79.31 ± 9.84 1.58 ± 0.39	-4.53 ± 4.01 4.76 ± 1.21	1.54 ± 1.03 0.33 ± 0.09	3.43 ± 2.43 0.42 ± 0.05
		BE GE .Plant.morphological.adaptations.. 9.50 ± 0.70 970.13 ± 125.37 5.31 ± 2.21 0.62 ± 0.26	54.56 ± 9.60 1.03 ± 0.38 TRY 19.21	-6.63 ± 2.13 2.84 ± 1.11 44	0.29 ± 0.40 0.42 ± 0.23	4.57 ± 1.62 0.51 ± 0.12
		PO SS root.metamorphoses. 8.21 ± 0.03 5.93 ± 4.73	780.40 ± 6.20 0.82 ± 0.13	83.15 ± 5.47 1.23 ± 0.47	-7.26 ± 1.62 2.97 ± 2.41	1.23 ± 0.38 0.44 ± 0.22	4.22 ± 1.11 0.46 ± 0.20
		GE CS .Root.carbon..C..content.per. 9.37 ± 0.05 789.33 ± 1.95 13.27 ± 6.64 0.79 ± 0.12	63.43 ± 8.16 2.00 ± 0.35 TRY 18.34	-5.91 ± 0.89 5.55 ± 2.09 42	0.36 ± 0.11 0.34 ± 0.07	5.22 ± 1.31 0.41 ± 0.08
		SS NO root.dry.mass. 7.54 ± 0.14 16.04 ± 6.39	734.67 ± 44.53 0.72 ± 0.13	55.89 ± 17.97 1.93 ± 0.50	-6.86 ± 0.57 6.78 ± 2.52	0.25 ± 0.79 0.37 ± 0.13	5.30 ± 0.96 0.44 ± 0.11
		CS Root_lignin_concentration 6.73 ± 0.05 Fun.rich	537.33 ± 13.43 Fun.even	56.64 ± 9.89 Fun.div GRooT 17.47	-5.34 ± 2.02 40	0.69 ± 0.49	3.61 ± 1.29
		NO IT .Root.carbon.nitrogen..C.N..ratio. 4.29 ± 1.56 821.11 ± 315.13 0.44 ± 0.13 0.75 ± 0.10	58.89 ± 13.90 0.022 ± 0.019 TRY 16.59	-6.48 ± 1.39 38	0.30 ± 0.49	3.68 ± 1.37
		Sand fraction 0.17 ± 0.12 Root_stele_diameter SW	pH 0.77 ± 0.17	mass(OS) 0.010 ± 0.009 GRooT 16.16	37
		(%) 24.94 ± 15.64 56.56 ± 5.63 66.66 ± 31.89 21.56 ± 12.27 42.58 ± 11.36 0.23 ± 0.19 0.28 ± 0.20 Root_stele_fraction IT SW FR BE PO FR BE PO 0.40 ± 0.16 Root_length_density_volume GE 0.24 ± 0.16 Root_P_concentration SS 0.20 ± 0.23 Specific_root_respiration CS 0.31 ± 0.16 .TRY 5.73 ± 0.78 50.18 ± 23.05 (g) 6.37 ± 0.85 129.74 ± 120.22 4.28 ± 0.22 47.16 ± 22.27 4.45 ± 0.42 53.97 ± 35.79 5.61 ± 0.69 92.95 ± 68.54 0.47 ± 0.27 0.014 ± 0.014 0.67 ± 0.18 GRooT 16.16 0.014 ± 0.014 0.59 ± 0.09 GRooT 15.72 0.010 ± 0.009 0.59 ± 0.26 GRooT 15.28 0.006 ± 0.003 0.71 ± 0.14 GRooT 13.97 0.005 ± 0.003 0.76 ± 0.10 0.004 ± 0.002 13.97	37 36 35 32 32
		GE SS NO .TRY 80.86 ± 1.87 4.24 ± 0.64 74.35 ± 39.84 30.45 ± 8.92 4.33 ± 0.30 67.66 ± 55.77 0.37 ± 0.10 0.63 ± 0.11 0.006 ± 0.003 13.54	31
		CS	43.73 ± 22.10	5.53 ± 0.21	36.06 ± 30.34
		NO Root_mass_density 60.89 ± 24.40	5.41 ± 0.37	25.60 ± 30.54 GRooT 13.10	30
		.Root.phosphorus..P..content.	TRY		12.66	29
		per.root.dry.mass.
		Net_nitrogen_uptake_rate		GRooT	12.23	28
		.TRY		12.23	28
		Rooting_depth		GRooT	10.48	24
		.Fine.root.carbon.nitrogen..C.N..ratio.	TRY		10.48	24
		.Fine.root.carbon..C..content.	TRY		9.61	22
		per.fine.root.dry.mass.
		Root_branching_density		GRooT	8.30	19
		Root_Ca_concentration		GRooT	8.30	19
		Root_K_concentration		GRooT	8.30	19
		.TRY		7.42	17
		Root_total_structural_carbohydrate_	GRooT	6.99	16
		Concentration
		Root_Mg_concentration		GRooT	6.11	14

Journal of Ecology: Confidential Review copy Journal of Ecology: Confidential Review copy

Journal of Ecology: Confidential Review copy

Research Fund (BOF) from Ghent University (01N02817). F.Z. received funding from the Swiss National Science Foundation (project 172198)

Acknowledgements

We thank Luc Willems and Greet De Bruyn for performing chemical analyses and Abdulwahhab Ghrairi for the texture analysis. We thank Evy Ampoorter, Haben Blondeel, Filip Ceunen, Kris Ceunen, Robbe De Beelde, Emiel De Lombaerde, Kent Hansson, Lionel Hertzog, Dries Landuyt, Pierre Lhoir, Sruthi M. Krishna Moorthy, Audrey Peiffer, Michael Perring, Matteo Tolosano, Sanne Van Den Berge, Lotte Van Nevel and Mia Vedel-Sørensen for providing help during the fieldwork campaign. We thank the handling editors and two anonymous reviewers for their constructive comments, which improved the manuscript significantly.

Funding: This work was supported by the Research Foundation Flanders (FWO) (K.D.P. ASP035-19, S.G. G0H1517N), the European Research Council (ERC) under the European Union's Horizon 2020 Research and Innovation Programme (ERC Starting Grant FORMICA 757833, 2018, http://www.formica.ugent.be) and the FWO Scientific research network FLEUR (http://www.fleur.ugent.be). T.V. received funding from the Special

Data availability

Data is available on Figshare. https://doi.org/10.6084/m9.figshare.c.5176136.v1

Authors' contributions PDF, SG, CM, PVG, KDP and KV conceived the ideas and designed methodology; all authors collected data; KDP analysed the data and led the writing of the manuscript, in close collaboration with PDF, PV, PS, SG, CM, and KV. All authors contributed critically to the drafts and gave final approval for publication.

Supplementary information

Supplementary information A: Soil texture and pH analysis To measure plot soil pH, five random subsamples of the mineral topsoil were taken (0-10 cm depth, diameter 30 mm) and pooled together. The samples were dried to constant weight at 40 °C for 48 h, ground and sieved over a 2 mm mesh. Then, they were analysed for pH-H 2 O by shaking a 1:5 ratio soil/H 2 O mixture for 5 min at 300 r.p.m. and measuring with an Orion 920A pH meter with a Ross sure-flow 8172 BNWP pH electrode model (Thermo Scientific Orion, USA). The same sampling procedure was followed for soil samples between 10-20 cm depth for texture analysis (% sand, silt and clay), which was performed by sieving and sedimentation with a Robinson-Köhn pipette according to ISO 11277 (2009).

Supplementary information B: GLMM analysis with study gradients from the experimental set-up as predictors (latitude, elevation, forest management type and distance to edge) Diversity metrics were modelled as a function of latitude, elevation, forest management type (dense, intermediate or open) and distance to forest edge, as these are the gradients constituting the experimental setup. All two-way interactions between these design variables were included in the models. Since the distance to the edge was designed in an exponential way, this predictor was log-transformed to meet model assumptions (log e ). Equation 2 represents the model structure.

Eq. 2

Supplementary Information C: GLMM analysis of floristic temperature of understorey communities To infer the temperature preference of species, we used the mean annual temperature of their distribution range. This value is based on a random sampling of 1000 mean annual temperature values across the distribution range of the species (De Frenne et al. 2013, Zellweger et al. 2020). We calculated the floristic temperature for the understorey communities, by taking the abundance weighted mean of the species' temperature values. For only seven out of 229 species, we were not able to calculate a temperature value (Agrostis mertensii, Helleborus viridis, Leopoldia comosa, Ranunculus bulbosus, Rubus plicatus, Taraxacum officinale, Valeriana excelsa).

In terms of data analyses, we applied the same modelling analysis which was performed on the diversity metrics (see Methods: data analyses) to the floristic temperature as response variable.

We found two environmental variables to drive the floristic temperature of the understorey communities: the mean annual temperature (MAT) and the winter minimum temperature offset. Unsurprisingly, for plots situated in regions with a higher mean annual temperature, the floristic temperature of the communities was significantly higher. More interestingly, the floristic temperature also increased when the microclimatic buffering of winter minimum temperatures was larger, i.e. when the minimum temperature was less cold inside the forest than in open terrain. We did not find a significant effect of microclimatic buffering of summer maximum temperatures or any other environmental variable.

Supplementary information D: Additional figures and tables

Tables