Simulating functional diversity of European natural forests along climatic gradients

We analyse how functional diversity (FD) varies across European natural forests to understand the effects of environmental and competitive filtering on plant trait distribution.


| INTRODUC TI ON
Functional diversity (FD) is a key control of the stability and adaptability of ecosystems under climate change (Yachi & Loreau, 1999).
Climatic conditions are the most important drivers of community assembly and generally constrain the relations between traits globally (Butler et al., 2017;Šímová et al., 2018). Gradients in trait expression, associated with climatic conditions, have been found at global scales (Díaz et al., 2016;Wright et al., 2017) and at the European scale by extrapolating site-specific plant trait data (Butler et al., 2017) and using forest inventories (Ruiz-Benito et al., 2017).
Recent studies investigated a mix of globally important physiological traits (e.g. SLA, WD, seed mass and leaf nitrogen content) and morphological traits (including maximum plant height and basal area) Ruiz-Benito et al., 2017), or separate the effects between those trait types (e.g. Madrigal-Gonzalez et al., 2016;Schneider et al., 2017). They advance our understanding of spatial pattern of plant traits at the landscape and continental scale.
Under given climate conditions, the distribution of multiple individual traits results from community assembly rules and opens a multi-dimensional trait space of functionally related traits (Chave et al., 2009;Díaz et al., 2016;Mason, Mouillot, Lee, & Wilson, 2005;Wright et al., 2004). To quantify how environmental and competitive filtering influence niche complementarity, multidimensional indices of FD are required to quantify occupation and overlap of niches . In addition, these indices should be scale independent and applicable at regional scale to investigate changes in FD along climatic gradients (Carmona, Bello, Mason, & Leps, 2016).
Many different indices have been used to describe the size of trait spaces, the distribution and the clustering of their trait combinations (e.g. Ratcliffe et al., 2016;Schneider et al., 2017). Villeger et al. (2008) suggested to describe these FD aspects using three independent indices: functional richness (FR), divergence (FDv) and evenness (FE).
Whereas FR quantifies the amount of occupied trait space, FDv and FE describe the distribution and abundance of trait combinations in a multidimensional trait space (Mason et al., 2005;Villeger et al., 2008).
Whereas FR describes the size of potentially available, functional space, in which niches can be occupied by plants, FDv quantifies the distribution of trait values, thus the degree of niche differentiation, likely the result from competitive exclusion (Garnier, Navas, & Grigulis, 2016;Mason et al., 2005). FE describes the regularity of trait distribution and points to resource-use efficiency within the occupied trait space. Lower FR could relate to lower capability of an ecosystem to buffer environmental stress, whereas lower FE and FDv could indicate reduced ecosystem resilience (Mason et al., 2005).
By dividing FD into richness, divergence and evenness, the mechanisms that link biodiversity to ecosystem functions and describe community assembly can be described (Mason, Bello, Mouillot, Pavoine, & Dray, 2013;Mason et al., 2005). Computing scale-independent FD indices from trait distributions (Carmona et al., 2016) allows investigating FD from community to meta-community scale and extrapolating them to the regional or continental scale. Huge efforts are under way to explore links among plant traits, vegetation composition and climate based on site data (Bruelheide et al., 2018;Díaz et al., 2016;Wright et al., 2017). They are complemented by dynamic global vegetation models (DGVM) with flexible or adaptive individual traits (e.g. Langan, Higgins, & Scheiter, 2017;Sakschewski et al., 2015) which explore and map the mechanisms between FD and ecosystem functions from local to regional scales. The interplay between flexible morphological and physiological traits within plant functional types (PFT, Lavorel et al., 2007;Prentice et al., 2007) in combination with the physiology and biogeochemistry of a DGVM allows analysing the effect of community assembly on ecosystem functions, for example, productivity and carbon storage, in forest ecosystems. Because these flexible-trait DGVMs vary plant traits for individual trees that belong to a specific PFT, inter-PFT as well as the intra-PFT trait diversity are captured which allows investigating effects of niche complementarity along climatic gradients.
The overall aim of this study is to investigate the interaction among climate, ecosystem functions and pattern of FD of forest ecosystems. Therefore, we adapted the DGVM Lund-Potsdam-Jena managed Land of flexible individual traits (LPJmL-FIT) to PFTs growing in strongly seasonal European climatic conditions, while the model previously has been successfully applied to tropical rainforests (Sakschewski et al., 2015). We break down the overall aim of the study into the following research questions: We focus on European natural forests ranging from broadleaved evergreen vegetation in the Mediterranean basin, to temperate forests and to boreal forests in northern Europe. Here, natural forests are defined as potential natural forests whose compositions and ecosystem functions (see Geller et al., 2017;Hooper et al., 2005 for definition) results from climate and soil conditions. Forests and other wooded land are usually defined following the percentage of woody cover (FAO, 2018). However, we denote forests hereafter as vegetation with a significant amount of biomass (>50 gC/m 2 ), at least 5% coverage of woody PFTs and a minimum mean tree height of 5 m.
To address these research questions, we check the validity of the adapted LPJmL-FIT model by assessing to what degree the model reconstructs observed (a) productivity and biomass, (b) trait distribution and (c) distribution of PFTs as a result of environmental and competitive filtering. We quantify for selected sites and on a Pan-European scale FR, FDv and FE of simulated physiological and morphological traits (cf. Schneider et al., 2017;Villeger et al., 2008). We expect FD to be relatively high in natural forests as was found in, for example, Schneider et al. (2017) and climatic stressors considered in LPJmL-FIT to be a strong environmental filter (Bernard-Verdier et al., 2012).
Our analysis of FD in natural forests can provide a reference state for restoring highly managed or degraded ecosystems and increase their diversity and stability in face of climate change (Mori, Lertzman, & Gustafsson, 2017). We focus here on the interaction of physiological and morphological plant traits with community assembly processes at the local and Pan-European scale to understand how FD emerges from those interactions.

| MATERIAL S AND ME THODS
We connected the leaf and stem economics approach as implemented in LPJmL-FIT (Sakschewski et al., 2015) with a phenology model (Forkel et al. (2014) to simulate potential natural vegetation in Europe under current climatic conditions. Herbaceous PFTs (C 3 and C 4 grasses) were simulated as in LPJmL . To allow for environmental and competitive filtering to take full effect within and across PFTs, we removed the bioclimatic limits that are used in most DGVMs to emulate biogeographic limitations of PFT occurrence Sitch et al., 2008). Phenology, leaf-economics (LES) and stem-economics (SES) traits (cf. Chave et al., 2009;Wright et al., 2004) were assigned to each individual tree sapling at establishment allowing any trait combination (everything-is-everywhere approach, see Figure 1 for trait-selection algorithm) whose competitiveness in a given climate then determines its survival and growth. The implemented LES traits include Specific Leaf Area (SLA, leaf area per unit leaf mass, mm 2 mg −1 ), Leaf Longevity (LL, average life span of leaves, in months), leaf nitrogen content (N leaf , leaf nitrogen content per leaf area, mg/g) and V cmax (maximum carboxylation rate of the RUBISCO enzyme per leaf area at 25°C, μmol CO 2 m −2 s −1 ). While the LES traits influence vegetation productivity, the SES trait Wood density (WD, wood dry mass per unit of green volume, g/cm 3 ) is related to biomass and tree mortality (for more details, see Sakschewski et al., 2015). Physiology, growth and mortality of trees within the forest patch were as described in Sakschewski et al. (2015). Applying this simulation framework to current European climate, environmental and competitive filtering result in site-specific trait distributions, productivity, biomass and in tree height.
LPJmL-FIT combines flexible individual traits with gap dynamics and plant physiology, hydrology and biogeochemistry. Being structured into vertical leaf layers every two meters, trees compete for light and water as they grow in size. The trait combination of each tree determines its competitive strength under given climate conditions at a given site, where several plant strategies can co-exist and form diverse communities in forest ecosystems. The suitability of the trait combinations to the local climate conditions determines trees' competitiveness, and thus which tree strategy is dominant from wet to dry tropical conditions (Sakschewski et al., 2015). In general, the model approach of LPJmL-FIT allows all tree strategies to establish everywhere at any time ("everything is everywhere"). Trees that underperform due to trait combinations less suitable for local climate and/or competitive conditions have a relatively high mortality probability (see Sakschewski et al., 2015 and 2.1 below). Under stable climate conditions, establishment and mortality lead to an equilibrium of tree abundance and community composition in which only those trees occur whose functional traits are suited to the local environmental conditions (climate and soil).
In the following, we describe the adjustments required to apply the LPJmL-FIT model to European natural forests.

| Adjusted LPJmL-FIT model
Lund-Potsdam-Jena managed Land of flexible individual traits was originally developed for Amazonian rainforests (Sakschewski et al., 2015). Adapting it to climatic conditions which are strongly seasonal and with a high intra-annual variability requires a number of adaptations and implementations for multiple, co-occurring PFTs to describe Mediterranean, temperate and boreal natural forests in Europe. The LES and SES approach, as implemented in LPJmL-FIT, has therefore been adapted to accommodate four tree PFTs, namely, "broad-leaved summer green" (BL-S), "broad-leaved evergreen" (BL-E), "Temperate needle-leaved evergreen tree" (T-NL) and "Boreal needle-leaved evergreen tree" (B-NL). Each PFT is based on an earlier implementation of these PFTs in LPJmL-4  which has been extensively evaluated . To account for the phenology of the Mediterranean, temperate and boreal forests under the influence of light, water and temperature stress, we implemented the phenology model of Forkel et al. (2014) into LPJmL-FIT. This model couples the phenological status of a tree, which ranges between 0 (complete senescence) and 1 (fully leaved), to the local climate. The actual value of the phenology status is determined by the product of four phenology functions, which depend on a set of PFT-specific parameters and the daily temperature, water stress and radiation. We calibrated the phenology parameter (Table S1) to yield a best possible PFT distribution that matches the spatial distribution of European natural vegetation from Bohn et al. (2007). In addition, leaf senescence now occurs immediately if the phenological status of a tree drops below 0.2, forcing a tree to rebuild the complete canopy leaf area in the next simulation year. Given the climate influence on leaf phenology changes along continental climate gradients, a continuous spectrum from summer to winter deciduousness for BL-S trees is captured in the model. The algorithm to combine the PFT parameter set of the LES-and SES-related traits (WD) with the phenology parameter is implemented as follows: When a new sapling establishes, the selection of the trait combination occurs in three steps. First, the PFT type is randomly chosen out of the four possible PFTs independent of its climate suitability or the present PFT composition ( Figure 1). Second, the selected PFT type defines the set of phenology parameters that are assigned to the tree sapling (see Table S1), and it defines the SLA range from which the SLA value is drawn from a uniform distribution. The SLA range differs for BL-E, T-NL, B-NL and BL-S (see Table S2). The functionally related plant traits leaf longevity (LL), leaf N and Vcmax (maximum carboxylation rate of Rubisco per leaf area) are assigned in a third step to the sapling following the LES approach as in Sakschewski et al. (2015), see anti-proportional to the actual tree cover and the mortality rate depends on carbon balance by the end of every year (see Schaphoff, von Bloh, et al., 2018 for details). Grasses are assigned phenology parameters (see Table S1), but are simulated without individual trait flexibility following the LPJmL4 modelling approach . Therefore, in this study model simulations will be used to evaluate trait distributions and FD regarding trees only.
Different from the LPJmL-FIT version of Sakschewski et al. (2015), the trade-off between SLA and LL has been adopted from LPJmL4 : where DM c denotes the dry matter carbon content of leaves DM c = 0.4763 and the parameter ∝= 2 ⋅ 10 −4 . The parameter 1 is set to 0.4 and 0 to 2.2 for broad-leaved PFTs, while 0 = 2.08 for both, T-NL and B-NL. All parameters in Equation (1) were obtained from Kattge et al. (2011). The parameters 1 and ∝ influence the steepness of the SLA-LL relation, whereas 1 alters the offset.
The mortality mort WD is coupled to its wood density WD by using the equation from King, Davies, Tan, and Noor (2006): and assigned to each tree individual at establishment (Sakschewski et al., 2015). Because no general mortality-WD relationship for tree species of temperate forests is currently available in the literature, we calibrated 1 and 2 to the locally observed biomass. Calibration was carried out at European sites still containing natural forests (Hainich National Park (NP) and Bialowieza NP) because the model simulates natural vegetation only. Following this calibration, we set the parameter 1 to −4.5 and 2 to −2.66 for the broad-leaved trees, and 1 = −2.66 and 2 = 0.255 for needle-leaved trees. The term mort WD is used as the maximum of the growth-efficiency mortality in LPJmL-FIT, meaning that trees with a low growth efficiency resulting from low productivity under unfavourable climate conditions have a higher mortality risk (Sakschewski et al., 2015).
In addition, we reduced the tree allometry parameter k ep to 1.5 (1.6 in  for needle-leaved trees to simulate realistic tree shapes and growth pattern of needle-in comparison to broad-leaved saplings: where k ep mediates between crown area (CA) and stem diameter (D).
The lower k ep , the CA for needle-leaved trees is reduced which then affects LAI: where C leaf is whole plant carbon investment in leaves (kg per tree) (cf. Sitch et al., 2003). Because of the lower SLA, needle-leaved trees then have to invest more leaf carbon in their first years to reach the same LAI compared with broad-leaved trees.
Fire is an important natural disturbance in European forest ecosystems (Naveh, 1990;Tinner et al., 1999). We applied the simple Glob-FIRM model (Thonicke, Venevsky, Sitch, & Cramer, 2001) as embedded in LPJmL4. Here, fire probability depends on soil moisture in the top soil layer and a fuel load threshold modelled by an exponential probability function calculated by the end of every year. We apply this probability to each patch separately. If a patch is burnt, every tree is ignited and survives with a PFT-specific fire resistance probability (cf. Thonicke et al., 2001) which is the same for all individuals that belong to the same PFT. In the patch burnt, all litter carbon is combusted completely.
In boreal forests, evergreen trees exhibit water stress in spring, when relatively high air temperature increases evaporative demand, while the soil is still frozen and thus limits the water availability in the soil. This water stress forces evergreen trees to shed their needles and making them less competitive against BL-S trees which might have their bud burst later in the year. To correctly balance competition between B-NL and BL-S in boreal forests, we increased the root-distribution factor root for B-NL trees from 0.943  to 0.965, which allows B-NL to reach deeper, non-frozen layers during spring thereby preventing them from leaf senescence.

| Model input data, simulation protocol and validation sites
To post-process and analyse our simulation data, we used R3.6.0 (R-Core-Team, 2019) and MATLAB 2019 (MATLAB, 2019). All R-packages used in this study including their citation are listed in Table S6.

| Simulation protocol
Model simulation starts from bare-ground and simulates a spin-up period of 500 years by recycling the first 30 years of the climate data set   Respective model output is then aggregated over all simulated patches within a grid cell.
To allow for a detailed analysis of plant-trait distribution, FD and productivity, we chose six different sites across Europe with near-natural forest stands and which cover a broad range of climates (Table S4). Site-specific simulations follow the same protocol as described above and were performed with 2,500 patches at each site to ensure a higher spatial coverage.

| Model validation
Simulated seasonal and intra-annual GPP is validated against observed and remotely sensed GPP data at six sites (Table S3) covering a climatic gradient (Table S4). We used monthly MODIS remote-sensing data (MOD17A2H) for the years 2004-2013 (Running, 2015) at six sites (see Table S1) and respective flux tower measurements from the Euroflux network for the Laegeren (CH-Lae, D'Odorico, 2014; Paul-Limoges, 2018) and Hainich NP (DE-Hai) (for general information, see Papale et al., 2006;Reichstein et al., 2005). Simulated maps of vegetation height and biomass were evaluated against remotely sensed products (Lefsky, 2007;Thurner et al., 2014). Details on the validation of Gross Primary Productivity (GPP), vegetation height and biomass are described in the Data S1.
Simulated plant trait distributions were compared against observed plant trait data from the TRY database (Kattge et al., 2011), see Data S1 for methods and data origin.

| Computation of FD indices
We quantified three complementary indices on multidimensional traits to describe FD, namely, FR, FDv and FE, following Villeger et al. (2008) and Schneider et al. (2017), where each point represents one tree individual (higher than 2 m) with its unique trait combination. In this study, this multidimensional trait space is based on SLA, LL, WD and tree height. While SLA, LL and WD influence productivity and biomass (Reich, 2014) and therefore point to competitive exclusion, tree height is regarded to describe niche differentiation (Garnier et al., 2016). We calculated the FD indices across and within PFTs to capture assembly processes across meta-communities.
Functional richness describes the extent of the occupied trait space and is calculated by the convex hull volume including all points in that trait space, which is normalized by the maximum possible trait volume. However, it implies that FR reacts strongly to outliers. FDv describes how far environmental niches are separated and indicate the intensity of competitive interactions, where FDv = 0 indicates convergent trait distribution due to strong environmental and competitive filtering (Mason et al., 2005;Villeger et al., 2008). To measure FDv in a multidimensional trait space, a sphere with radius dG centred in the trait cloud is calculated. FDv then quantifies how points (trait combinations of trees) scatter relative to the surface of the sphere (see eqs. 5-7 in Data S1 and Figure S1). If all points are located on the sphere, FDv becomes unity independent on the respective distribution on the sphere.
The more the index decreases, the wider the points are spread around (inside and outside) the sphere. FDv therefore quantifies how the occupied niches are separated. FE describes how regularly points are distributed in the trait space, that is, how efficient available resources are used through niche occupation. FE is based on the minimum spanning tree (MST) linking all points in the trait space in such a way that the sum of all branches becomes minimal (see eqs. 11-13 in Data S1 and Figure S1). Therefore, FE increases when a) the points are evenly distributed, that is, having equal branch length; or b) the trait combinations of the trees are equidistant in the trait space (Villeger et al., 2008). Further details on FE, FDv and FR quantification are provided in Data S1.
Each index was computed by considering all trees as evenly weighted. Before calculation, traits were normalized to their minimum and maximum observable values in the TRY sites in Europe, thus ranging between 0 and 1. For each grid cell we checked that no dimensional reduction was required by using the function "dbFD" of the R-package "FD" performing a principal coordinates analysis (Laliberté et al., 2014). Due to constrained computation capacity, we calculated all FD indices separately in groups of 50 patches in each grid cell (for which 1,000 patches were simulated in total) and aggregated them to the grid level by using the arithmetic mean. To visualize the stochastic uncertainty of the model, we calculated the coefficient of variation (COV) of each index in a grid cell out of the groups of 50 patches (n = 20). As all of these groups in a cell received the same climate data, the COV can be seen as a measure for the stochastic uncertainty of the model.

| Climate influence on trait distribution, productivity and tree height
Environmental and competitive filtering allows those trees to establish and survive whose trait combinations are suitable for local climate conditions in the forest patches simulated by LPJmL-FIT.
Unsuitable or less suitable trait combinations lead to a low growth efficiency and are therefore outcompeted. The combination of climate suitability and competiveness has the effect that the continuous, uniform distribution with which the model is initialized results in a normal trait distribution at the local scale. The resulting trait distributions therefore emerge from the LPJmL-FIT modelling framework ( Figure 1). We compare simulated SLA and WD against TRY observations for BL-S, BL-E trees and needle-leaved evergreen, that is, B-NL and T-NL trees ( Figure 2). Simulated mean trait distributions match the TRY observation reasonably well for both, SLA and WD, for BL-S, BL-E and the two needle-leaved PFTs (dashed lines in Figure 2). Simulated ranges of SLA, however, are smaller than the original trait range (see Table S2) and smaller than observed SLA ( Figure 2). LPJmL-FIT simulates a mean SLA of 16.11 with a standard deviation of ±1.25 mm 2 /mg for BL-S compared to 14.95 ± 6.09 documented in TRY. The simulated range for SLA is also smaller for BL-E (LPJmL-FIT: 9.31 ± 1.13 mm 2 /mg; TRY: 6.95 ± 2.72 mm 2 /mg) and the needle-leaved evergreen (LPJmL-FIT: 6.11 ± 0.95 mm 2 /mg; TRY: 6.95 ± 2.72 mm 2 /mg). Simulated ranges for WD are quite close to observed ranges in TRY for BL-S, slightly smaller for BL-E, but broader for the needle-leaved PFTs (Figure 2, bottom row).  Table S3, Figure S2 and Data S1 for details on the evaluation methods and results, sites are described in Table S4). The simulated competition between tree individuals results in closed forest cover and corresponding high biomass storage in temperate and boreal forests ( Figure S3). We found simulated and reduced natural forests to few remaining small areas (Ellis et al., 2013), which further complicates the evaluation of simulated potential natural vegetation.

| Climate influence on trait distributions within and across PFTs
Spatial distribution of simulated fractional cover of each PFT result from the PFT-specific phenology and functional trait combinations, which determine the suitability of each parameter set to the climate and soil conditions in a given grid cell (cf. Figure 1).  Table S2). BL-E cover a similar temperature range as BL-S, but occur across a wider SLA range between  Figure 3a; Figure S5 shows PFT-specific SLA maps, Figure S8 maps the COV of SLA).  Figure S6 for PFT-specific WD maps and Figure S9 for maps displaying the COV of WD).

| FD emerging from climate and plant competition
The calculation of all diversity indices is based on the initial four-dimensional trait space out of SLA, LL, WD and tree height. However, for visualization, we remapped the trait space from four dimensions to a three-dimensional trait space composed of SLA, WD and tree height. We plot the position of each tree individual in the trait space for the climatologically different sites Seitseminen, Laegern and Dundo (Figure 4, left column). The occupied trait space forms the hypervolume, that is, the FR (shown in grey-blue in Figure 4, 2nd column) for each site. The sphere around the centre of gravity (grey surface and green cross, respectively, shown in Figure 4, 3rd column) illustrates site-specific FDv, whereas FE is quantified from the MST ( Figure 4, last column). Table S4 shows the site-specific FD indices for the six sites.
The wide bi-modal distribution of SLA between BL-S and B-NL trees in Seitseminen increases the trait space, that is, FR, whereas in Laegern and Dundo simulated SLA distributions show narrower bi-modal distributions or even converge (see density distribution in Figure 4, left column). Niche separation (FDv) and regularity of niche occupation (FE) are more comparable across the three sites (Table S4).
FDv is highest in Dundo because points in trait space lay closer to the surface of the sphere compared with Seitseminen and Laegern (notice points outside the sphere in Seitseminen and Laegern). Compared with Laegern, we find slightly higher FDv in Seitseminen because of the divergent SLA distribution. Niche occupation is less regularly distributed (FE) in Dundo compared with Seitseminen and Laegern because the trait space of B-NL trees is less occupied in Dundo, leading to larger path length in between points of this PFT (Data S1, Figure S1). This leads to more irregular distances in between points, which lowers FE.

F I G U R E 4
Distribution of simulated trees in a three-dimensional trait space (Height, WD and SLA). We plot density (far left) and functional diversity indices for 2,500 simulated patches (FR centre left; FDv centre right and FE (far right) for three sites (Seitseminen (top row), Lägern (middle row) and Dundo (bottom row)). Ranges of the SLA and WD axes correspond to the maximum trait range across all woody PFTs used in the simulation (see Table S2). Seitseminen, Laegern and Dundo represent boreal, mixed temperate and Mediterraneantype forests respectively. Each dot represents a trait combination of one tree larger than 5 m while the colour indicates its PFT type. Plant Functional Types are as followed: Broad-leaved summer green (BL-S), Broad-leaved Evergreen (BL-E), Boreal needle-leaved evergreen tree (B-NL) and Temperate needle-leaved evergreen tree (T-NL). FE, functional evenness; FR, functional richness; PFT, plant functional type; SLA, specific leaf area; WD, wood density When calculating the FD indices for each PFT trait space separately and for each site, FDv and FE are similar to the corresponding across-PFT values (Table S5 and Data S1 for computation of within-PFT FD indices). Within each PFT, diverse functional strategies co-exist, that is, the intensity of plant competition and regularity of niche occupation are comparable to the one across PFTs. Specifically, FDv is high at Seitseminen for B-NL trees compared with the overall FDv because the point cloud is clearly separated along the tree-height niche axis and shows a wider WD distribution compared with the other sites ( Figure 4, left column). However, intra-PFT FR is three orders of magnitude lower than FR across PFTs. This is mainly caused by a much lower realized trait space for LL (e.g. BL-S) and SLA (e.g. B-NL, see Figure 4), whereas the intra-PFT range of tree heights and WDs is similar to that between PFTs. In summary, environmental and competitive filtering influence niche occupation in a similar way within as well as across PFTs.

| FD at the European scale
At the European scale, spatial gradients in FD indices are relatively small and a few spatially distinct patterns stand out ( Figure 5).  Figure S7).
Functional divergence is high in natural forests reaching values between 0.68 and 0.82 (Figure 5b). Where needle-leaved and broadleaved trees co-exist (cf. Figure S4), FDv is higher, that is, in boreal and mountain forests, and in southern Mediterranean forests.
Where only one PFT dominates, FDv is lower ( Figure S7).

| Effect of environmental and competitive filtering on trait distribution and productivity
The adapted LPJmL-FIT model is capable of reproducing observed GPP with a small modelling error and high correlation with observed data (see Data S1). Biomass and plant height follow the spatial distribution of previous publications Healey et al., 2015), although the comparison is limited by the long- Simulated ranges for SLA (BL-S and BL-E) and WD (BL-E and needle-leaved trees) are smaller than observed (Figure 2 and Table   S2). Heat and drought stress as well as light availability (seasonal, vertical canopy structure) influence growing conditions of trees.  Figure S3). We have included fire disturbance in LPJmL-FIT, which simulates high fire activity in the Mediterranean forests, less in boreal and low fire activity in temperate forests, but other disturbance agents such as windthrow or frost could also further diversify plant strategies, and thus increase the width of the simulated trait distribution. Implementing these disturbances into LPJmL-FIT would be expected to decrease the competitiveness of the BL-S trees in boreal and mountainous forests, and thus further restrict their spatial extent.
Specific model functions determine the model outcome. The phenology parameter (Table S1)  With more SLA and LL measurements available for temperate forests, this relationship could be adjusted and contribute to further reduce modelling error of the simulated trait distribution (Figure 2).

| Community-assembly effects on traits and FD across Europe
Community assembly at a site results from dispersal or migration, Simulated SLA cover a wide temperature and precipitation range, separated by SLA ranges as observed for the four PFTs in Europe.
Where cold temperatures and light increasingly limit productivity, that is, in the boreal zone, simulated SLA for BL-S increases, meaning trees with extremely high SLA and short LL survive cold winters and grow in short summers. The shift towards higher WD in southern Europe indicates a local adaptation to seasonally dry Mediterraneantype climate. In general, the increasingly dry climate filters higher WD for BL-E and T-NL which indicates better adaptation to drought.  (Garnier et al., 2016;Mason et al., 2005;Villeger et al., 2008).  Figure S3). We find that such bi-modal changes in FDv reflect shifts between competitive exclusion, linked to physiological traits (SLA, WD, and LL), and niche differentiation, linked to morphological traits (tree height).
Functional evenness quantifies the regularity of the distribution in the trait space, that is, niche occupation (cf. Mason et al., 2005;Villeger et al., 2008) and could be interpreted as a lower utilization of resources due to the more irregular occupation of the trait space, that is, environmental niches (cf. Mason et al., 2005). However, there are few studies which have investigated the changes in FE along climatic gradients. Recent studies focussed on changes in FE along disturbance gradients (e.g. Mouillot, Graham, Villeger, Mason, & Bellwood, 2013) and at the local scale (Pakeman, 2011

S U PP O RTI N G I N FO R M ATI O N
Additional supporting information may be found online in the Supporting Information section.
How to cite this article: Thonicke