This is the accepted version of the journal article: Thonicke, Kirsten; Billing, Maik; Von Bloh, Werner; [et al.]. «Simulating func- tional diversity of European natural forests along climatic gradients». Journal of biogeography, Vol. 47, Issue 5 (May 2020), p. 1069-1085. DOI 10.1111/jbi.13809 This version is available at https://ddd.uab.cat/record/257068 under the terms of the license 1 Simulating functional diversity of 2 European natural forests along climatic 3 gradients 4 to be submitted to Journal of Biogeography 5 Running title: Functional diversity of European forests 6 7 Authors: Kirsten Thonicke1,*, Maik Billing1,2,*, Werner von Bloh1, Boris Sakschewski1, Ülo 8 Niinemets3, Josep Peñuelas4, Hans Cornelissen5, Yusuke Onoda6, Peter van Bodegom7, 9 Michael Schaepman8, Fabian Schneider9, Ariane Walz2 10 1Research Domain 1 “Earth System Analysis”, Potsdam Institute for Climate Impact 11 Research (PIK), Member of the Leibniz Association, P.O. Box 60 12 03, D-14412 Potsdam, 12 Germany 13 2 Institute of Environmental Science and Geography, University of Potsdam, Karl-Liebknecht- 14 Str. 24-25, 14476 Potsdam-Golm, Germany 15 3 Estonian University of Life Sciences, Kreutzwaldi 1, 51006 Tartu, Estonia 16 4 Global Ecology Unit CREAF-CSIC-UAB, Center for Ecological Research and Forestry 17 Applications (CREAF) - National Research Council (CSIC) Edifici C, Universitat Autonoma 18 de Barcelona, 08193 Bellaterra 19 5 SystemsEcology, Department of Ecological Science, Vrije Universiteit, De Boelelaan 1085, 20 1081 HV Amsterdam, the Netherlands 21 6 Divisionof Environmental Science and Technology, Graduate School of Agriculture, Kyoto 22 University, Oiwake, Kitashirakawa, Kyoto, 606-8502, Japan 23 7 Institute of Environmental Sciences, Department Environmental Biology, Leiden University, 24 Einsteinweg 2, 2333 CC Leiden 25 8 Remote Sensing Laboratories, Dept. of Geography, University of Zurich, 26 Winterthurerstrasse 190, CH–8057 Zurich, Switzerland 27 9Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, 28 Pasadena, CA 91011, USA 29 * Kirsten Thonicke and Maik Billing should be considered joint first author 30 Abstract 31 1 32 Aim 33 We analyse how functional diversity varies across European natural forests to understand 34 the effects of environmental and competitive filtering on plant trait distribution. 35 Location 36 Forest ecosystems in Europe from 11°W to 36°E and 29.5°N to 62°N. 37 Taxon 38 Pinaceae, Fagaceae and Betulaceae, Oleaceae, Tiliaceae, Aceraceae, Leguminosae 39 (unspecific). 40 Methods 41 We adopted the existing Dynamic Global Vegetation Model Lund-Potsdam-Jena managed 42 Land of flexible individual traits (LPJmL-FIT) for Europe by eliminating both bioclimatic limits 43 of plant functional types (PFTs) and replacing prescribed values of functional traits for PFTs 44 with emergent values under influence of environmental filtering and competition. We 45 quantified functional richness (FR), functional divergence (FDv) and functional evenness 46 (FE) in representative selected sites and at Pan-European scale resulting from simulated 47 functional and structural trait combinations of individual trees. While FR quantifies the 48 amount of occupied trait space, FDv and FE describe the distribution and abundance of trait 49 combinations, respectively, in a multidimensional trait space. 50 Results 51 LPJmL-FIT reproduces spatial PFTs and local trait distributions and agrees well with 52 observed productivity, biomass and tree height of European natural forests. The observed 53 site-specific trait distributions and spatial gradients of traits of the leaf- and stem-resource 54 economics spectra coincide with environmental filtering and the competition for light and 55 water in environments with strong abiotic stress. Where deciduous and needle-leaved trees 56 co-occur, e.g. in boreal and mountainous forests, the potential niche space is wide (high 57 FR), and extreme ends in the niche space are occupied (high FDv). We find high FDv in 58 Mediterranean forests where drought increasingly limits tree growth, thus niche 59 differentiation becomes more important. FDv decreases in temperate forests where a cold 60 climate increasingly limits growth efficiency of broadleaved summergreen trees, thus 61 reducing the importance of competitive exclusion. Highest FE was simulated in wet Atlantic 62 and southern Europe which indicated relatively even niche occupation and thus high 63 resource-use efficiency. 64 Main conclusions 65 We find functional diversity resulting from both environmental and competitive filtering. Pan- 66 European FR, FDv and FE demonstrate the influence of climate gradients and intra- and 67 inter-PFT competition. The indices underline a generally high functional diversity of natural 68 forests in Europe. Co-existence of functionally diverse trees across PFTs emerges from 69 alternative (life-history) strategies, disturbance and tree demography. 2 70 1. INTRODUCTION 71 Functional diversity (FD) is a key control of the stability and adaptability of ecosystems under 72 climate change (Yachi and Loreau, 1999). Abiotic conditions (e.g. climate, soil) as well as 73 biotic processes (e.g. competition) determine plant community assembly (Kunstler et al., 74 2016; Ratcliffe et al., 2017; Ruiz-Benito et al., 2017), and thus the functional diversity of 75 communities (Lavorel & Garnier, 2002a; Cadotte et al., 2011; Naeem et al., 2012). 76 Ecosystem FD has been measured based on the diversity of morphological and/or 77 physiological plant traits (Villeger et al., 2008), which are linked to plant productivity, 78 transpiration or nutrient cycling, and thus to ecosystem functions (cf. Lavorel & Garnier, 79 2002b). For instance, specific leaf area (SLA) and leaf nitrogen content have been linked to 80 plant productivity, and stem traits such as wood density (WD) to carbon storage and 81 evapotranspiration (Funk et al., 2017) and, among others, determine species tolerance of 82 environmental stresses (Hallik et al., 2009) and control competitive interactions between 83 individual plants (Kunstler et al., 2016). 84 Under given climate conditions, the distribution of multiple individual traits results from 85 community assembly rules and opens a multi-dimensional trait space of functionally related 86 traits (Wright et al., 2004; Mason et al., 2005; Chave et al., 2009; Díaz et al., 2016). To 87 quantify how environmental and competitive filtering influence niche complementarity, 88 multidimensional indices of functional diversity are required to quantify occupation and 89 overlap of niches (Mason & de Bello, 2013). In addition, these indices should be scale- 90 independent and applicable at regional scale to investigate changes in FD along climatic 91 gradients (Carmona et al., 2016). 92 Climatic conditions are the most important drivers of community assembly and generally 93 constrain the relations between traits globally (Butler et al., 2017; Šímová et al., 2018). 94 Strong gradients in trait expression, associated to climatic conditions, have been found at 95 global scales (Díaz et al., 2016; Wright et al., 2017) and at the European scale by 96 extrapolating site-specific plant trait-data (Butler et al. 2017) and using forest inventories 97 (Ruiz-Benito et al. 2017). Recent studies investigate a mix of globally important physiological 98 traits (e.g., SLA, WD, seed mass and leaf nitrogen content) and morphological traits 99 (including maximum plant height and basal area) (Ratcliffe et al., 2016; Ruiz-Benito et al., 100 2017), or separate the effects between those trait types (e.g. Madrigal-Gonzalez et al., 2016; 101 Schneider et al., 2017). They advance our understanding of spatial pattern of plant traits at 102 the landscape and continental scale. 103 Many different indices have been used to describe the size of trait spaces, the distribution 104 and the clustering of their trait combinations (e.g. Ratcliffe et al., 2016; Schneider et al., 105 2017). Villéger et al. (2008) suggested to describe these FD aspects using three 106 independent indices: functional richness (FR), divergence (FDv) and evenness (FE). While 107 FR quantifies the amount of occupied trait space, FDv and FE describe the distribution and 108 abundance of trait combinations in a multidimensional trait space (Mason et al., 2005; 109 Villeger et al., 2008). While FR describes the size of potentially available, functional space, 110 in which niches can be occupied by plants, FDv quantifies the distribution of trait values, thus 3 111 the degree of niche differentiation, likely the result from competitive exclusion (Mason et al., 112 2005; Garnier et al., 2016). FE describes the regularity of trait distribution and points to 113 resource-use efficiency within the occupied trait space. Lower FR could relate to lower 114 capability of an ecosystem to buffer environmental stress, whereas lower FE and FDv could 115 indicate reduced ecosystem resilience (Mason et al., 2005). 116 By dividing functional diversity into richness, divergence and evenness, the mechanisms that 117 link biodiversity to ecosystem function and describe community assembly can be described 118 (Mason et al., 2005; Mason et al., 2013). Computing scale-independent FD indices from trait 119 distributions (Carmona et al., 2016) allows investigating FD from community to meta- 120 community scale and extrapolating them to the regional or continental scale. Huge efforts 121 are under way to explore links between plant traits, vegetation composition and climate 122 based on site data (Díaz et al., 2016; Wright et al., 2017; Bruelheide et al., 2018). They are 123 complemented by dynamic global vegetation models (DGVM) with flexible or adaptive 124 individual traits (e.g., Sakschewski et al., 2015; Langan et al., 2017) which explore and map 125 the mechanisms between functional diversity and ecosystem function from local to regional 126 scales. The interplay between flexible morphological and physiological traits within plant 127 functional types (PFT, Lavorel et al., 2007; Prentice et al., 2007) in combination with the 128 physiology and biogeochemistry of a DGVM allows analysing the effect of community 129 assembly on ecosystem function, e.g. productivity and carbon storage, in forest ecosystems. 130 Because these flexible-trait DGVMs vary plant traits for individual trees that belong to a 131 specific PFT, inter-PFT as well as the intra-PFT trait diversity are captured which allows 132 investigating effects of niche complementarity along climatic gradients. 133 The overall aim of this study is to investigate the interaction between climate, ecosystem 134 function and pattern of FD of forest ecosystems. We focus on European natural forests 135 ranging from broadleaved evergreen vegetation in the Mediterranean basin, to temperate 136 forests and to boreal forests in northern Europe. Here, natural forests are defined as 137 potential natural forests whose compositions and ecosystem functions (see Mitchell, 2002; 138 Hooper et al., 2005; Geller et al., 2017 for definition) results from climate and soil conditions. 139 Forests and other wooded land are usually defined following the percentage of woody cover 140 (FAO, 2018). However, we denote forests hereafter as vegetation with a significant amount 141 of biomass (>50 gC/m²), at least 5% coverage of woody PFTs and a minimum mean tree 142 height of 2m. For this study we adapted the DGVM LPJmL-FIT to PFTs growing in strongly 143 seasonal European climatic conditions, while the model previously has been successfully 144 applied to tropical rainforests (Sakschewski et al., 2015). We break down the overall aim of 145 the study into the following research questions: 146 1. What is the role of environmental and competitive filtering on trait distributions and 147 productivity of Mediterranean, temperate and boreal natural forests? 148 2. How does functional diversity emerge from climate and plant competition at the local 149 and pan-European scale? 150 3. How does functional diversity vary between and across European natural forests? 151 In order to address these research questions, we check the validity of the adapted LPJmL- 152 FIT model by assessing to what degree the model reconstructs observed (a) productivity and 4 153 biomass, (b) trait distribution and (c) distribution of PFTs as a result of environmental and 154 competitive filtering. We quantify for selected sites and on a Pan-European scale FR, FDv 155 and FE of simulated physiological and morphological traits (cf. Villeger et al., 2008; 156 Schneider et al., 2017). We expect functional diversity to be generally high in natural forests 157 and climatic stressors considered in LPJmL-FIT to be a strong environmental filter. 158 Our analysis of functional diversity in natural forests can provide a reference state for 159 restoring highly managed or degraded ecosystems and increase their diversity and stability 160 in face of climate change (Liang et al., 2016; Mori et al., 2017). We focus here on the 161 interaction of physiological and morphological plant traits with community assembly 162 processes at the local and Pan-European scale to understand how functional diversity 163 emerges from those interactions. 164 2. MATERIALS AND METHODS 165 We connect the leaf and stem economics approach as implemented in LPJmL-FIT 166 (Sakschewski et al., 2015) with a phenology model (Forkel et al. (2014) to simulate potential 167 natural vegetation in Europe under current climatic conditions. Herbaceous PFTs (C3 and C4 168 grasses) are simulated as in LPJmL (Schaphoff et al., 2018b). To allow for environmental 169 and competitive filtering to take full effect within and across PFTs, we removed the 170 bioclimatic limits that are used in most DGVMs to emulate biogeographic limitations of PFT 171 occurrence (Sitch et al., 2008; Schaphoff et al., 2018b). Phenology, leaf-economics (LES) 172 and stem-economics (SES) traits (cf. Wright et al., 2004; Chave et al., 2009) are assigned to 173 each individual tree sapling at establishment allowing any trait combination (everything-is- 174 everywhere-approach, see Fig. 1 for trait-selection algorithm) whose competitiveness in a 175 given climate then determines its survival and growth. Physiology, growth and mortality of 176 trees within the forest patch are as described in Sakschewski et al. (2015). Applying this 177 simulation framework to current European climate, environmental and competitive filtering 178 result in site-specific trait distributions, productivity, biomass and in tree height. 179 2.1 Adjusted LPJmL-FIT model 180 LPJmL-FIT combines flexible individual traits with gap dynamics and plant physiology, 181 hydrology and biogeochemistry. Being structured into vertical leaf layers every two meters, 182 trees compete for light and water as they grow in size. The trait combination of each tree 183 determines its competitive strength under given climate conditions at a given site, where 184 several plant strategies can co-exist and form diverse communities in forest ecosystems. 185 The competitive interactions between individual trees through the suitability of their trait 186 combinations to the given climate, spatial distribution of traits change along respective 187 climatic gradients (Sakschewski et al., 2015). In general, the model approach of LPJmL-FIT 188 allows all tree strategies to establish everywhere at any time (“everything is everywhere”). 189 However, in equilibrium only those trees survive that are best adapted to the local 190 environmental conditions (climate and soil). 5 191 LPJmL-FIT was originally developed for Amazonian rainforests (Sakschewski et al., 2015). 192 Adapting it to climatic conditions which are strongly seasonal and with a high intra-annual 193 variability requires a number of adaptations and implementations for multiple, co-occurring 194 PFTs to describe Mediterranean, temperate and boreal natural forests in Europe. The LES 195 and SES approach, as implemented in LPJmL-FIT, has therefore been adapted to 196 accommodate four tree PFTs, namely “broad-leaved summergreen” (BL-S), “broad-leaved 197 evergreen” (BL-E), “temperate needle-leaved” (T-NL) and “boreal needle-leaved trees” (B- 198 NL). Each PFT is based on an earlier implementation of these PFTs in LPJmL-4 (Schaphoff 199 et al., 2018b) which has been extensively evaluated (Schaphoff et al., 2018a). To account 200 for the phenology of the Mediterranean, temperate and boreal forests under the influence of 201 light, water and temperature stress, we implemented the phenology model of Forkel et al. 202 (2014) into LPJmL-FIT. This model couples the phenological status of a tree, which ranges 203 between zero (complete senescence) and one (fully leafed), to the local climate. The actual 204 value of the phenology status is determined by the product of four phenology functions, 205 which depend on a set of PFT-specific parameters and the daily temperature, water stress 206 and radiation. We calibrated the phenology parameter (Table S1) to yield a best possible 207 PFT distribution that matches the spatial distribution of European natural vegetation from 208 Bohn et. al (2007). Additionally, leaf senescence now occurs immediately if the phenological 209 status of a tree drops below 0.2, forcing a tree to rebuild the complete canopy leaf area in 210 the next simulation year. Given the climate influence on leaf phenology, BL-S trees imply a 211 continuous spectrum to winter-deciduous trees. 212 The algorithm to combine the PFT parameter set of the LES and SES-related traits with the 213 phenology parameter is implemented as follows: When a new sapling is established, the 214 selection of the trait combination occurs in three steps. First, the PFT type is randomly 215 chosen out of the four possible PFT independent of its climate suitability or the present PFT 216 composition (Fig. 1). Second, the selected PFT type defines the set of phenology 217 parameters that are assigned to the tree sapling (see Suppl. Table S1), and it defines the 218 SLA range from which the SLA value is drawn from a uniform distribution. The SLA range 219 differs for BL-E, T-NL, B-NL and BL-S (see Table S2). The functionally related plant traits 220 leaf longevity (LL), leaf N and Vcmax (maximum carboxylation rate of Rubisco per leaf area) 221 are assigned in a third step to the sapling following the LES approach as in Sakschewski et 222 al. (2015), see Fig. 1. Both, SLA and WD are drawn from a continuous uniform distribution 223 for each individual tree for which the PFT-specific range is derived from the TRY v4.0 224 database (see https://www.try-db.org/ Kattge et al., 2011) considering only sites located 225 inside our study area. Hence, the parameter set of each newly established sapling in the 226 forest patch contains trait values drawn from the LES and SES as well as the phenology 227 parameters. Grasses are integrated in the model as homogeneous layers using the 228 remaining radiation at the bottom of every forest patch for photosynthesis. C3 and C4 grasses 229 compete with trees for water only. The establishment rate of grass is anti-proportional to the 230 actual tree cover and the mortality rate depends on carbon balance by the end of every year 231 (see Schaphoff et al. 2018 for details). Grasses are assigned phenology parameters (see 232 Table S1), but are simulated without individual trait flexibility following the LPJmL4 modelling 6 233 approach (Schaphoff et al., 2018b). Therefore, in this study model simulations will be used to 234 evaluate trait distributions and functional diversity regarding trees only. 235 Different to the LPJmL-FIT version of Sakschewski et al. (2015), the trade-off between SLA 236 and LL has been adopted from LPJmL4 (Schaphoff et al., 2018b): 𝑆𝐿𝐴 𝛽0 −log10 ( ∙𝐷𝑀𝑐 ) ∝ 237 𝐿𝐿 = 10 𝛽1 (1) 238 where 𝐷𝑀𝑐 denotes the dry matter carbon content of leaves (𝐷𝑀𝑐 = 0.4763) and the 239 parameter ∝ = 2 ∙ 10−4 . The parameter 𝛽1 is set to 0.4 and 𝛽0 to 2.2 for broadleaved PFTs, 240 while 𝛽0 = 2.08 for both, T-NL and B-NL. All parameters in equation (1) were obtained from 241 Kattge et al. (2011). The parameters 𝛽1 and ∝ influence the steepness of the SLA-LL 242 relation, whereas 𝛽1 alters the offset. 243 The mortality mortWD is coupled to its wood density WD by using the equation from King et 244 al. (2006): 245 𝑚𝑜𝑟𝑡𝑊𝐷 = 10𝛼1+𝛼2/𝑊𝐷 (2) 246 and assigned to each tree individual at establishment (Sakschewski et al., 2015). Because 247 no general mortality-WD relationship for tree species of temperate forests is currently 248 available in the literature, we calibrated 𝛼1 and 𝛼2 to the locally observed biomass. 249 Calibration was carried out at European sites still containing natural forests (Hainich National 250 Park (NP) and Bialowieza NP), because the model simulates natural vegetation only. 251 Following this calibration, we set the parameter 𝛼1 to -4.5 and 𝛼2 to -2.66 for the broadleaved 252 trees, and 𝛼1 = −2.66 and 𝛼2 = 0.255 for needle-leaved trees. The term mortWD is used as 253 the maximum of the growth-efficiency mortality in LPJmL-FIT, meaning that trees with a low 254 growth efficiency resulting from low productivity under unfavourable climate conditions have 255 a higher mortality risk (Sakschewski et al., 2015). 256 Additionally, we reduced the tree allometry parameter kep to 1.5 (1.6 in Schaphoff et al., 257 2018b) for needle-leaved trees to simulate realistic tree shapes and growth pattern of 258 needle- in comparison to broadleaved saplings: 259 𝐶𝐴 ~ 𝐷 𝑘𝑒𝑝 (3) 260 where kep mediates between crown area (CA) and stem diameter (D). The lower kep, the CA 261 for needle-leaved trees is reduced which then affects LAI: 𝐶𝑙𝑒𝑎𝑓 ∙𝑆𝐿𝐴 262 𝐿𝐴𝐼 = (4), 𝐶𝐴 263 where 𝐶𝑙𝑒𝑎𝑓 is whole plant carbon investment in leaves (kg per tree) (cf. Sitch et al., 2003). 264 Because of the lower SLA, needle-leaved trees then have to invest more leaf carbon in their 265 first years to reach the same LAI compared to broadleaved trees. 266 Fire is an important natural disturbance in European forest ecosystems (Naveh, 1990; 267 Tinner, 1999). We applied the simple Glob-FIRM model (Thonicke et al., 2001) as 268 embedded in LPJmL4. Here, fire probability depends on soil moisture in the top soil layer 7 269 and a fuel load threshold modelled by an exponential probability function calculated by the 270 end of every year. We apply this probability to each patch separately. If a patch is burnt, 271 every tree is ignited and survives with a PFT-specific fire resistance probability (cf. Thonicke 272 et al., 2001) which is the same for all individuals that belong to the same PFT. In the patch 273 burnt, all litter carbon is combusted completely. 274 In boreal forests, evergreen trees exhibit water stress in spring, when relatively high air 275 temperature increases evaporative demand, while the soil is still frozen and thus limits the 276 water availability in the soil. This water stress forces evergreen trees to shed their needles 277 and making them less competitive against BL-S trees which might have their bud burst later 278 in the year. To correctly balance competition between B-NL and BL-S in boreal forests, we 279 increased the root-distribution factor 𝛽𝑟𝑜𝑜𝑡 for B-NL trees from 0.943 (Schaphoff et al., 280 2018b) to 0.965, which allows B-NL to reach deeper, non-frozen layers during spring thereby 281 preventing them from leaf senescence. 282 2.2 Model input data, simulation protocol and validation sites 283 Climate and soil data input 284 LPJmL-FIT uses air temperature [°C], precipitation [mm/d] and radiation (short-wave down 285 and long-wave net radiation [W/m²]) of the combined dataset of the WATCH (Weedon, 2011) 286 and WFDEI (Weedon, 2014) datasets at daily resolution on a 0.5° x 0.5° longitude-latitudinal 287 grid. This climate data set is based on the reanalysis of ERA-Interim, where precipitation 288 was bias-corrected using the Global Precipitation Climatology Centre data set (GPCC, 289 Schneider et al., 2011). The climate data range from 1901 to 2013 with WFD covering 1901 290 to 1978 and WFDEI-GPCC is used from 1979 onwards. The atmospheric CO2 concentration 291 is held constant at 296 ppm over the whole time period. 292 Soil texture is needed as model input and was taken from the Harmonized World Soil 293 Database (HWSD) version 1.2 (Nachtergaele & Montanarella, 2009). The soil depth was 294 kept constant at 2m for all grid cells. Our simulation domain covers Europe from 11°E to 295 36°W and 29.5°N to 62°N. 296 Simulation protocol 297 Model simulation starts from bare-ground and simulates a spin-up period of 500 years by 298 recycling the first 30 years of the climate data set (1901-1930) to bring natural vegetation 299 composition (here individual trees with their individual trait combinations) and all living and 300 dead carbon fluxes into equilibrium with the spin-up climate. We then performed a transient 301 run simulating potential natural vegetation until the end of 2013, i.e. without land use. For the 302 European simulation domain, 1000 forest patches being equivalent to 10 ha of forest area 303 are simulated in each grid cell where all patches receive the same climate data and the 304 same soil data as model input. Respective model output is then aggregated over all 305 simulated patches within a grid cell. 8 306 To allow for a detailed analysis of plant trait distribution, functional diversity and productivity, 307 we chose six different sites across Europe with near-natural forest stands and which cover a 308 broad range of climates (Table S4). Site-specific simulations follow the same protocol as 309 described above and were performed with 2500 patches at each site to ensure a higher 310 spatial coverage. 311 Model validation 312 Simulated seasonal and intra-annual GPP is validated against observed and remotely 313 sensed GPP data at six sites (Table S3) covering a climatic gradient (Table S4). We used 314 monthly MODIS remote sensing data (MOD17A2H) for the years 2004-2013 (Running, 315 2015) at six sites (see Table S1) and respective flux tower measurements from the Euroflux 316 network for the Laegeren (CH-Lae, D'Odorico, 2014; Paul-Limoges, 2018) and Hainich NP 317 (DE-Hai) (for general information, see Reichstein et al., 2005; Papale et al., 2006). 318 Simulated maps of vegetation height and biomass were evaluated against remotely sensed 319 products (Lefsky, 2007; Thurner, 2014). Details on the validation of GPP, vegetation height 320 and biomass are described in the SI. Simulated plant trait distributions were compared 321 against observed plant trait data from the TRY data base (Kattge et al., 2011), see SI for 322 methods. 323 2.3 Computation of functional diversity indices 324 We quantified three complementary indices on multi-dimensional traits to describe FD, 325 namely functional richness (FR), functional divergence (FDv), and functional evenness (FE) 326 following Villéger et al. (2008) and Schneider et al. (2017), where each point represents one 327 tree individual (higher than 2m) with its unique trait combination. In this study, this 328 multidimensional trait space is based on SLA, LL, WD and tree height. While SLA, LL and 329 WD influence productivity and biomass (Reich, 2014) and therefore point to competitive 330 exclusion, tree height is regarded to describe niche differentiation (Garnier et al., 2016). We 331 calculated the FD indices across and within PFTs to capture assembly processes across 332 meta-communities. 333 FR describes the extent of the occupied trait space and is calculated by the convex hull 334 volume including all points in that trait space, which is normalized by the maximum possible 335 trait volume. However, it implies that FR reacts strongly to outliers. FDv describes how far 336 environmental niches are separated and indicate the intensity of competitive interactions, 337 where FDv = 0 indicates convergent trait distribution due to strong environmental and 338 competitive filtering (Mason et al., 2005; Villeger et al., 2008). To measure FDv in a 339 multidimensional trait space, a sphere with radius ̅̅̅̅ 𝑑𝐺 centred in the trait cloud is calculated. 340 FDv then quantifies how points (trait combinations of trees) scatter relative to the surface of 341 the sphere (see eqs. 5-7 in SI and Fig. S1). If all points are located on the sphere, FDv 342 becomes unity independent on the respective distribution on the sphere. The more the index 343 decreases, the wider the points are spread around (inside and outside) the sphere. FDv 344 therefore quantifies how the occupied niches are separated. FE describes how regularly 345 points are distributed in the trait space, i.e. how efficient available resources are used 9 346 through niche occupation. FE is based on the minimum spanning tree (MST) linking all 347 points in the trait space in such a way, that the sum of all branches becomes minimal (see 348 equ. 11-13 in SI and Fig. S1). Therefore, FE increases when a) the points are evenly 349 distributed, i.e. having equal branch length; or b) the trait combinations of the trees are 350 equidistant in the trait space (Villeger et al., 2008). Further details on FE, FDv and FR 351 quantification are provided in the SI. 352 Each index was computed by considering all trees as evenly weighted. Before calculation, 353 traits were normalized to their minimum and maximum observable values in the TRY sites in 354 Europe, thus ranging between 0 and 1. For each grid cell we checked that no dimensional 355 reduction was required by using the function “dbFD” of the R-package “FD” performing a 356 principal coordinates analysis (Laliberté, 2014). Due to constrained computation capacity, 357 we calculated all FD indices separately in groups of 50 patches in each grid cell (for which 358 1000 patches were simulated in total) and aggregated them to the grid level by using the 359 arithmetic mean. To visualize the stochastic uncertainty of the model, we calculated the 360 coefficient of variation (COV) of each index in a grid cell out of the groups of 50 patches 361 (n=20). Since all of these groups in a cell received the same climate data, the COV can be 362 seen as a measure for the stochastic uncertainty of the model. 363 3. RESULTS 364 3.1 Climate influence on trait distribution, productivity and tree 365 height 366 Environmental and competitive filtering allows those trees to establish and survive whose 367 trait combinations are suitable for local climate conditions in the forest patches simulated by 368 LPJmL-FIT. Unsuitable or less suitable trait combinations lead to a low growth efficiency and 369 are therefore outcompeted. The combination of climate suitability and competiveness has 370 the effect that the continuous, uniform distribution with which the model is initialized results 371 in a normal trait distribution at the local scale. The resulting trait distributions therefore 372 emerge from the LPJmL-FIT modelling framework (Fig.1). We compare simulated SLA and 373 WD against TRY observations for BL-S, BL-E trees and needle-leaved evergreen, i.e. B-NL 374 and T-NL trees (Fig. 2). Simulated mean trait distributions match the TRY observation 375 reasonably well for both, SLA and WD, for BL-S, BL-E and the two needle-leaved PFTs 376 (dashed lines in Fig. 2). Simulated ranges of SLA, however, are smaller than the original trait 377 range (see Table S2) and smaller than observed SLA (Fig.2). LPJmL-FIT simulates a mean 378 SLA of 16.11 with a standard deviation of +-1.25 mm²/mg for BL-S compared to 14.95 +- 379 6.09 documented in TRY. The simulated range for SLA is also smaller for BL-E (LPJmL-FIT: 380 9.31 +- 1.13 mm²/mg; TRY: 6.95 +- 2.72 mm²/mg) and the needle-leaved evergreen 381 (LPJmL-FIT: 6.11 +- 0.95 mm²/mg; TRY: 6.95 +- 2.72 mm²/mg). Simulated ranges for WD 382 are quite close to observed ranges in TRY for BL-S, slightly smaller for BL-E, but broader for 383 the needle-leaved PFTs (Fig. 2, bottom row). 10 384 The combined effects of environmental filtering and plant competition for light and water in 385 LPJmL-FIT also result in reasonable seasonal and interannual productivity (gross primary 386 productivity, GPP) as observed on six selected sites. The model quality is shown by a high 387 Pearson’s R (≥0.89) and low NMSE (0.03 to 0.29; see Table S3, Fig. S2 and SI for details on 388 the evaluation methods and results, sites are described in Table S4). The simulated 389 competition between tree individuals results in closed forest cover and corresponding high 390 biomass storage in temperate and boreal forests (Fig. S3). We found simulated vegetation 391 height and biomass to compare well against remotely sensed observations and local in-situ 392 data, although the comparison of natural forest and actual vegetation is limited as remote- 393 sensing products detect properties of actual vegetation cover which are influenced by 394 current land-use and forest management (see SI for details on the evaluation method and 395 results). Centuries of land clearing, agricultural and forestry have greatly changed land cover 396 and reduced natural forests to few remaining small areas (Ellis et al., 2013), which further 397 complicates the evaluation of simulated potential natural vegetation. 398 3.2 Climate influence on trait distributions within and across 399 PFTs 400 Spatial distribution of simulated fractional cover of each PFT result from the PFT-specific 401 phenology and functional trait combinations, which determine the suitability of each 402 parameter set to the climate and soil conditions in a given grid cell (cf. Fig.1). Tree 403 individuals with trait combinations adapted to current climate are the most competitive and 404 most productive in sites with limited environmental stress severity (while their productivity 405 may decrease at stressed sites (Zhang et al., 2018), and thus cover larger proportions of a 406 given grid cell (Fig. S4). The most suitable combination, to shed leaves under cold and/or 407 dry climate conditions, results in BL-S dominating central Europe, even though it also occurs 408 - albeit at much smaller fractions- in the boreal and the Mediterranean forests. B-NL 409 dominates northern Europe, and T-NL the Mediterranean basin, where it co-occurs with BL- 410 E. Several PFTs co-occur in the Mediterranean forests (3 tree PFTs and C3 grasses) while 411 temperate forests in lowland Europe are dominated by just one tree PFT (Fig. S4). Note that 412 in all cases, the tree individuals still vary in their trait combination within each PFT. 413 Abiotic conditions, here aggregated to MAT [°C] and MAP [mm], are strongly linked to tree 414 establishment (Fig. 3). The climate space occupied by trees across all PFTs converges 415 towards higher MAT and lower MAP. Intra-PFT variation in SLA values decreases with 416 warmer and drier conditions. Most SLA variation, however, happens between PFTs that 417 occupy different parts of the climate space (Fig. 3a). The SLA of BL-S varies from values 418 around 13 mm²/mg under warmer and drier climate conditions to >20 mm²/mg in colder 419 climate conditions (<7°C) with a wide range in precipitation (500 to >2000 mm MAP). BL-S 420 co-exists with B-NL in cold/wet climate conditions, with BL-E in warmer and increasingly 421 drier climate conditions (>10°C MAT and 300 to 2000 mm MAP). Although many trait 422 combinations are possible under the “everything-is-everywhere” approach of LPJmL-FIT, 423 BL-S with SLA values lower than 13 mm²/mg do not occur despite a possible minimum of 7 424 mm²/mg (see Table S2). BL-E cover a similar temperature range as BL-S, but occur across 11 425 a wider SLA range between a MAT of 8 to 17°C, although the possible SLA range of BL-E is 426 among the smallest of all simulated PFTs. Again, this realized trait space separates these 427 PFTs clearly from the two needle-leaved PFTs, T-NL and B-NL, where T-NL increasingly 428 dominates at MAT>10°C and MAP<1000mm with SLA ranging from 5 to 8 mm²/mg, whereas 429 B-NL shows lowest SLA values at a MAT below 10°C (see Fig. 3a; Fig. S5 shows PFT- 430 specific SLA maps). 431 Simulated mean WD clearly separates along climate gradients across all PFTs (Fig. 3b). 432 Low WDs are simulated in cold climate conditions, and WD increases with increasing MAT 433 and decreasing MAP. B-NL and BL-S cover the WD space between 0.5 and <0.7 g/cm³ 434 below 10°C MAT. WD around 0.7 g/m³ are found in cold (<10°C MAT) and drier climate 435 conditions (<1000 mm MAP), with similar WD values found in warmer (10-12°C MAT) and 436 wetter (>1000 mm MAP) climate conditions. Highest WD values (>1.1 g/cm³) are simulated 437 for T-NL, BL-S and BL-E with >15°C MAT and <1000 mm MAP (Fig. 3b). Higher WD allows 438 slow plant growth and lowers tree mortality risk (eq. 2 in 2.1) which explains why high WD 439 are simulated across all PFTs under dry climate conditions (see also Fig. S6 for PFT-specific 440 WD maps). 441 3.3 Functional diversity emerging from climate and plant 442 competition 443 The calculation of all diversity indices is based on the initial 4-dimensional trait space out of 444 SLA, LL, WD and tree height. However, for visualisation, we remapped the trait space from 445 four dimensions to a three dimensional trait space composed of SLA, WD and tree height. 446 We plot the position of each tree individual in the trait space for the climatological different 447 sites Seitseminen, Laegern and Dundo (Fig. 4, left column). The occupied trait space forms 448 the hypervolume, i.e. the FR (shown in grey-blue in Fig. 4, 2nd column) for each site. The 449 sphere around the centre of gravity (grey surface and green cross, respectively shown in 450 Fig. 4, 3rd column) illustrates site-specific FDv, while FE is quantified from the minimum 451 spanning tree (Fig. 4, last column). Table S4 shows the site-specific FD indices for the six 452 sites. 453 The wide bi-modal distribution of SLA between BL-S and B-NL trees in Seitseminen 454 increases the trait space, i.e. FR, whereas in Laegern and Dundo simulated SLA 455 distributions show narrower bi-modal distributions or even converge (see density distribution 456 in Fig. 4, left column). Niche separation (FDv) and regularity of niche occupation (FE) are 457 more comparable across the three sites (Table S4). FDv is highest in Dundo, because points 458 in trait space lay closer to the surface of the sphere compared to Seitseminen and Laegern 459 (notice points outside the sphere in Seitseminen and Laegern). Compared to Laegern, we 460 find slightly higher FDv in Seitseminen, because of the divergent SLA distribution. Niche 461 occupation is less regularly distributed (FE) in Dundo compared to Seitseminen and 462 Laegern, because the trait space of B-NL trees is less occupied in Dundo, leading to larger 463 path length in between points of this PFT (SI, Fig. S1). This leads to more irregular distances 464 in between points, which lowers FE. 12 465 When calculating the FD indices for each PFT trait space separately and for each site, FDv 466 and FE are similar to the corresponding across-PFT values (Table S5 and SI for 467 computation of within-PFT FD indices). Within each PFT, diverse functional strategies co- 468 exist, i.e. the intensity of plant competition and regularity of niche occupation is comparable 469 to the one across PFTs. Specifically, FDv is high at Seitseminen for B-NL trees compared to 470 the overall FDv, because the point cloud is clearly separated along the tree-height niche axis 471 and shows a wider WD distribution compared to the other sites (Fig. 4, left column). 472 However, intra-PFT FR is three orders of magnitude lower than FR across PFTs. This is 473 mainly caused by a much lower realized trait space for LL (e.g. BL-S) and SLA (e.g. B-NL, 474 see Fig. 4), while the intra-PFT range of tree heights and WDs is similar to that between 475 PFTs. In summary, environmental and competitive filtering influence niche occupation in a 476 similar way within as well as across PFTs. 477 3.4 Functional diversity at the European scale 478 At the European scale, spatial gradients in FD indices are relatively small and a few spatially 479 distinct patterns stand out (Fig. 5). FR increases with the number of PFTs present through 480 which the size of the trait space increases (Fig. 5a, compare Fig. S4). Higher FR is found in 481 mountain areas throughout the continent, in boreal forests but also on the British Isles. 482 Where B-NL and BL-S occupy distant parts of the trait space, FR reaches its maximal values 483 of 0.03 to 0.04, very much alike in the Seitseminen site (compare Fig.4, top row). In contrast, 484 lowest FR values are computed for areas where one PFT is dominant, especially in the 485 lowland areas of temperate forests which are dominated by BL-S trees (Fig. 5a). The 486 coefficient of variation (COV) for FR is high in temperate and alpine forests, where mean FR 487 is low, because variability increases where the mean of a variable is close to zero. On the 488 contrary, low values are found in Mediterranean and boreal forests (Fig. S7). 489 Functional Divergence is high in natural forests reaching values between 0.68 and 0.82 (Fig. 490 5b). Where needle-leaved and broad-leaved trees co-exist (cf. Fig. S4), FDv is higher, i.e. in 491 boreal and mountain forests, and in southern Mediterranean forests. Where only one PFT 492 dominates, FDv is lower (0.7-0.73 compared to >0.75), e.g. in lowland temperate and 493 Mediterranean forests. In the transition zone to boreal forests or mountain forests, trait 494 distributions of BL-S-dominated forests further converges (FDv ~0.68). Here, intensity of 495 filtering increases (FDv converges) when climate conditions reduce growth efficiency of BL-S 496 trees. Further north, growing conditions for B-NL are more suitable, allowing establishment 497 of another plant strategy causing an increase in FDv (Fig. 5b). In contrast, in the southern 498 Mediterranean forests, needle-leaved and broad-leaved trees are smaller (Fig. S3) and in 499 competition with grasses for water, thus FDv is higher again (~0.76, see Fig. 5b). COV of 500 FDv is generally very low with maximum values found in some mountainous and boreal 501 transitional areas (Fig. S7). 502 Functional Evenness follows a different spatial pattern across Europe than the other two 503 indices (Fig. 5c). Lower FE values (~0.72) are found in lowland temperate forest in central 504 and eastern Europe as well as in the transition zone to the boreal forest in north-eastern 505 Europe, where BL-S dominate but still grow in competition with B-NL trees (Fig. S3). FE 13 506 increases where only one NL-PFT dominates (e.g. British Isles, Norway) and is high in wet 507 Atlantic and in southern Europe, where high environmental stress results in efficient 508 resource use. Here, COV of FE corresponds to the spatial gradient of mean FE, maximum 509 COV is found where mean FE is lower (Fig. S7). 510 4. DISCUSSION 511 4.1 Effect of environmental and competitive filtering on trait 512 distribution and productivity 513 The adapted LPJmL-FIT model is capable of reproducing observed GPP with a small 514 modelling error and high correlation with observed data (see SI). Biomass and plant height 515 follow spatial distribution of previous publications (Thurner, 2014; Healey et al., 2015), 516 although the comparison is limited by the long-term land-use history in Europe which 517 restricts the comparability between simulated and observed data and may explain the 518 discrepancy found. Simulated SLA and WD are in close agreement with observed TRY data. 519 These results show that the new version of LPJmL-FIT reproduces the spatial PFT and local 520 trait distributions as well as the productivity and biomass of European natural forests. Even 521 though a preference for measuring short-lived, broad-leaved species can influence 522 measured trait data which could overestimate, e.g., SLA (Sandel et al., 2015), such bias is 523 found to be small for European trees and should have little influence on our results because 524 LPJmL-FIT clearly separates between broad-leaved and needle-leaved trees and focuses on 525 trees only. The model is capable of simulating potential natural vegetation without 526 prescribing the spatial extent of PFTs via bioclimatic limits and without prescribing the 527 functional traits in question (e.g. SLA, LL, WD). Therefore, it enables to investigate the 528 interaction between environmental and competitive filtering in European natural forests in an 529 unprecedented manner. 530 Simulated ranges for SLA (BL-S and BL-E) and WD (BL-E and needle-leaved trees) are 531 smaller than observed (Fig. 2 and Table S2). Heat and drought stress as well as light 532 availability (seasonal, vertical canopy structure) influence growing conditions of trees. These 533 climate factors act as environmental filter favouring plant strategies that are adapted to the 534 local climate. Tree growth efficiency influences tree competiveness and its capacity to reach 535 the top-layer in the canopy and gain most light in temperate and boreal forests. Therefore, 536 competitive filtering is the second driver of trees surviving and growing in a forest patch. 537 Additional processes not yet explicitly captured by LPJmL-FIT might explain the remaining 538 trait variability and include: i) dispersal and adaptive responses (incl. phenotypic plasticity of 539 traits), ii) nutrient availability (e.g. nitrogen limitation), iii) variable rooting strategies and 540 disturbances other than fire (Douma et al., 2012; Van Bodegom et al., 2012), and iv) trade- 541 offs between different trait combinations and species capacities to tolerate multiple 542 environmental stresses simultaneously (Niinemets & Valladares, 2006; Laanisto & 543 Niinemets, 2015). While the influence of adaptive responses on trait distribution is perhaps 544 difficult to measure, including nutrient availability would introduce another niche axis and 14 545 perhaps differentiate the trade-off with SLA in the LES better. Variable rooting strategies 546 could alleviate water stress in seasonally dry environments and increase growth efficiency, 547 hence tree height of Mediterranean forests (cf. Fig. S3). We have included fire disturbance in 548 LPJmL-FIT, which simulates high fire activity in the Mediterranean forests, less in boreal and 549 low fire activity in temperate forests, but other disturbance agents such as wind-throw or 550 frost could also further diversify plant strategies, and thus increase the width of the simulated 551 trait distribution. Implementing these disturbances into LPJmL-FIT would be expected to 552 decrease the competitiveness of the BL-S trees in boreal and mountainous forests, and thus 553 further restricting their spatial extent. 554 Specific model functions determine the model outcome. The phenology parameter (Table 555 S1) influence GPP, growth efficiency, and thus the performance of a tree in the forest patch. 556 In exchange with climate conditions, the phenology parameters determine the spatial 557 distribution of the PFTs, including grasses. The trade-off between SLA and LL defines the 558 carbon cost for leaves (or needles) under specific climate conditions which then determine 559 the length of the growing season, and thus vegetation productivity. With more SLA and LL 560 measurements available for temperate forests this relationship could be adjusted and 561 contribute to further reduce modelling error of the simulated trait distribution (Fig.2). The 562 WD-dependent mortality function (see equ. 2) influences simulated biomass, more data on 563 WD and tree mortality would help to improve parameterisation of this function. Further 564 reductions in the model error of simulated WD distribution, tree height and biomass could be 565 expected. 566 Plant-trait validation would profit from more and better resolved plant trait measurements in 567 natural forests in Europe. In this study, we had to aggregate several TRY sites and also to 568 merge T-NL and B-NL to validate simulated trait distribution. Ideally, trait measurements at 569 FLUXNET sites, where we can evaluate both productivity and plant traits, would allow high 570 resolution and in-depth model evaluation. With these consistently measured data available, 571 model evaluation could be extended to water fluxes and biomass as well as stand structure 572 and their related plant traits. Alternatively, remotely sensed traits which are already available 573 at the Laegern site (cf. Schneider et al., 2017) or emerging from global remote-sensing 574 mission (cf. Jetz et al., 2016) could be used in the future. However, cross-validation 575 experiments are required to compare spatially continuous traits observed from remote 576 sensing to traits simulated by DGVMs using flexible-individual traits, but approaches are 577 underway (cf. Garonna, 2018; Lausch, 2018). Simulated trait maps for SLA and WD can be 578 validated against interpolated observed trait maps (Butler et al., 2017; Šímová et al., 2018), 579 even though observed trait distributions include those from highly managed forests. 580 Interpolated trait maps, which combine plant trait data with actual species’ presence data, 581 can lead to a potential bias in data-model comparison because species are planted in 582 managed forests. Generating trait maps which use potential natural species distribution 583 would improve such data-model comparison because they would be closer to the spatial trait 584 distribution that would result from environmental and competitive filtering. Model-data 585 comparison could be improved with a) more trait data being extracted in natural vegetation 586 and extrapolated using species distribution data to the continental scale; and b) allowing for 587 a better overlap of observed vs. simulated plant traits such as leaf N (observation available 15 588 but simulated map missing) or wood density (observation missing, but simulated map 589 available). 590 4.2 Community-assembly effects on traits and functional 591 diversity across Europe 592 Community assembly at a site results from dispersal or migration, environmental filtering of 593 trait combinations adapted to local climate and finally plant competition (Bernard-Verdier et 594 al., 2012). In the adapted LPJmL-FIT model, trait combinations are drawn from observed 595 ranges for SLA and WD. It principally means that any conceivable trait combination can 596 establish everywhere at any time for every PFT without explicitly considering trait inheritance 597 or seed dispersal. Currently, surviving trait combinations in the model represent those trees 598 whose traits are best suited to local climate and the competitive conditions in the established 599 tree community. Therefore, environmental filtering interacts with competitive filtering and 600 reproduces observed trait distributions and productivity of natural forests. Scaling up to the 601 whole study region, SLA and WD gradients emerge along climatic gradients and reflect 602 variability in trait ranges. 603 Simulated SLA cover a wide temperature and precipitation range, separated by SLA ranges 604 as observed for the 4 PFTs in Europe. Where cold temperatures and light increasingly limit 605 productivity, i.e. in the boreal zone, simulated SLA for BL-S increases, meaning trees with 606 extremely high SLA and short LL survive cold winters and grow in short summers. The shift 607 towards higher WD in southern Europe indicates a local adaptation to seasonally dry 608 Mediterranean-type climate. Generally, the increasingly dry climate filters higher WD for BL- 609 E and T-NL which indicates better adaptation to drought. 610 The identified spatial FD patterns reflect the combined effect of environmental and 611 competitive filtering. The influence of climate on surviving plant strategies and their 612 coexistence are reflected in the occupied trait space. FR is high where trait clusters of 613 different PFTs are distant and thus cover a large volume of the trait space, which is the case 614 in high-elevation areas and the boreal zone, where BL-S and B-NL co-exist. Here, BL-S 615 adapt to high SLA values and lower LL, because of shorter vegetation periods leading to a 616 wider functional gap in these two traits between BL-S and B-NL trees. FR is small where 617 only BL-S trees dominate the patches of temperate forests (Fig. 5). Here, trait variation is 618 small in each of the traits used, specifically LL, thus resulting in a small trait space, whereas 619 the ranges of WD and tree height do not change substantially. In the Mediterranean region, 620 where the trait space covered by each PFT is larger, FR is comparable to temperate forests. 621 Although the coexistence of up to 4 PFTs would lead to a relatively high FR in leaf traits in S- 622 Europe, lower tree heights limit trait space and counteract the effect of diverse leaf 623 strategies. In the southern Mediterranean, where tree growth reaches its limits, FR 624 decreases because few niches can be occupied under the dry climate conditions. 625 Functional divergence describes the diversity of co-existing functional strategies. A high FDv 626 illustrates a high degree of niche differentiation, thus high competitive exclusion or intense 627 competitive interaction (Mason et al., 2005; Villeger et al., 2008; Garnier et al., 2016). The 16 628 larger the influence of environmental filtering through abiotic stress, the more the trait 629 distribution converges (Bernard-Verdier et al., 2012). Maps of FDv (Fig.5) therefore illustrate 630 the relative strength of environmental vs. competitive filtering on trait composition. Where 631 climate allows needle- and broad-leaved trees to co-exist (in mountains, boreal forests, wet 632 temperate and Mediterranean forests), FDv is high because of the distant SLA distribution of 633 BL-S and B-NL. In Mediterranean and temperate climates, SLA distributions converge more, 634 and FDv declines slightly to values around 0.7 (Fig.5b). Interestingly, FDv declines further in 635 the transition from temperate forests to boreal and to mountain forests, respectively. Being 636 dominated by BL-S, the performance of these trees declines as climate gets colder, thus the 637 abiotic stress increases which leads to a less divergent distribution. Further north (or at 638 higher elevation) the growing conditions for B-NL improve, thus their dominance increases 639 and increases FDv again. Physiological traits (SLA, WD, LL) converge at the transition from 640 temperate to boreal forests, indicating competitive exclusion. In contrast, in southern 641 Mediterranean forests, increasing drought stress reduces tree height and favours a wide 642 SLA range and higher WD (Fig. 3, S5). In the northern mountain ranges, T-NL, BL-S and B- 643 NL co-exist, each with different SLA and WD. T-NL and B-LE with SLA and WD close to 644 each other dominate the lowland and coastal area, resulting in lower FDv. The distributions 645 diverge again further south with increasing drought stress, promoting co-occurring 646 alternative strategies (Fig. 5b, S3). We find that such bi-modal changes in FDv reflect shifts 647 between competitive exclusion, linked to physiological traits (SLA, WD, LL), and niche 648 differentiation, linked to morphological traits (tree height). 649 Functional evenness quantifies the regularity of the distribution in the trait space, i.e. niche 650 occupation (cf. Mason et al., 2005; Villeger et al., 2008) and could be interpreted as a lower 651 utilization of resources due to the more irregular occupation of the trait space, i.e. 652 environmental niches (cf. Mason et al., 2005). However, there are few studies which have 653 investigated the changes of FE along climatic gradients. Recent studies focussed on 654 changes in FE along disturbance gradients ( e.g., Mouillot et al., 2013) and at the local scale 655 (Pakeman, 2011). We find high FE in areas where dry and cold climate conditions limit tree 656 growth and under wet (Atlantic) conditions. Increasing disturbances (fire) and climatic stress 657 (drought) increase resource-use efficiency as suggested by Pakeman (2011). In highly 658 stressed environments, no strategy can dominate, thus functional traits are more evenly 659 distributed. Low FE is found in mixed temperate forests and in the transition from temperate 660 to boreal forests, where BL-S and B-NL occupy distant parts of the trait space. At the Pan- 661 European scale, FE does not follow a unique climatic gradient which could be interpreted as 662 habitat filtering as suggested by Pakeman (2011). In our study, it is rather a combination of 663 environmental and competitive filtering within and across PFTs. 664 We computed FR, FDv and FE from four traits following the approach of Villeger et al. (2008) 665 as implemented for the Laegern site in Schneider et al. (2017). The FD indices change with 666 the traits considered and the dimensionality of the trait space. At the Laegern site, we 667 compute FR values one order of magnitude lower than published in Schneider et al. (2017) 668 because we calculate FR for a 4- instead of a 3-dimensional trait space. In contrast, we 669 derived comparable FDv and slightly lower FE values for the Laegern site (Schneider et al., 670 2017). A direct comparison between remotely sensed and simulated FD based on 3 traits is 17 671 of interest and relies on directly comparable trait spaces and indices. Quantification of FD 672 indices strongly depend on the type of traits (physiological vs. morphological) and the 673 context of their interpretation, e.g. niche differentiation or ecosystem function to which the 674 considered traits relate. Comparing site-specific FD indices based on similar traits where 675 also trait observations are available would be a great step forward to evaluate site-specific 676 FD and its changes along climatic gradients at the regional as well as continental scale. 677 Model-data comparison can help to further investigate the importance of spatial scale or 678 abiotic gradients at the landscape scale, where DGVMs with flexible individual traits such as 679 LPJmL-FIT can be used to test and explore respective FD hypotheses to advance 680 biodiversity-ecosystem function theory (Hisano et al., 2018). 681 Trees occupy trait ranges within each PFT allowing different PFTs to co-exist as a result of 682 functional diversity, tree demography and disturbance in all biomes. Site-specific water 683 availability and temperature have been shown to influence functional composition in 684 temperate forests emphasizing the importance to investigate functional diversity along 685 climate gradients (Zhang et al., 2018). 686 We find high functional diversity in European natural forests not only for temperate forests as 687 found by Liebergesell et al. (2016) but also for boreal and Mediterranean forests. The 688 simulated community assembly of European biomes can be explained by the interplay 689 between environmental and competitive filtering: I) In the boreal zone and in high-mountain 690 areas (e.g. Carpathian Mountains, Alps or the Pyrenees) the growing season is short, and 691 needle-leaved trees (low SLA) which keep their leaves longer (high LL) but grow relatively 692 fast due to low WD, are the dominant tree growth strategy mostly as a result of 693 environmental filtering. B-NL trees have high resistance to cold temperatures and their 694 evergreen strategy allows them to fully exploit the short growing season (see Table 1). In 695 contrast, broadleaved trees (BL-S), which are more prone to cold temperatures and have to 696 unfold their thin, short-lived leaves at the start of the growing season, have a disadvantage 697 under this climate. Nevertheless over a wide range within the boreal zone, B-NL trees are 698 accompanied by BL-S trees. This means, in this range environmental and competitive 699 filtering are not strong enough to select for B-NL trees only and both growing strategies can 700 co-exist. Competitive filtering though pushes BL-S towards lower LL and higher WD as in the 701 temperate zone, which minimizes the overlap of trait values between the two PFTs (Fig. 1). 702 Due to the BL-S and B-NL co-existence a wide trait space (high FR) is covered in which trait 703 combinations are mostly evenly spaced (high FE) and divergent (high FDv), allowing for 704 relatively high resource-use efficiency. The same effect applies further South towards the 705 temperate zone, where B-NL trees can still be found in mixed temperate forests. 706 II) In the temperate zone, climate supports a longer growing season, higher productivity and 707 biomass. Here, BL-S trees dominate due to competitive filtering. Their tolerance of warmer 708 temperatures together with low carbon investment in short-lived and thin leaves (Table 1) 709 makes this PFT very competitive along a large climate gradient. Moreover, towards warmer 710 and drier climate conditions, the WD of BL-S trees increases in order to survive this climatic 711 stress keeping this PFT dominant. The clear PFT dominance decreases FR, and leads to 712 intermediate FE and FDv. Towards the Mediterranean zone new tree growth strategies (BL- 18 713 E, T-NL) enter the forest community with milder winter and drier summer changing 714 environmental filtering effects. 715 III) In the Mediterranean zone forests are stressed by drought and warm temperatures, 716 therefore limitations of growing conditions vary and environmental filtering plays an important 717 role. BL-E and T-NL trees with low SLA and high WD are adapted to seasonally dry 718 summers and can therefore co-exist in complementary and regularly spaced niches which 719 confirms previous findings (Carnicer et al., 2013). Whereas BL-E trees tolerate warm 720 temperatures and drought, T-NL trees are specialists for warm, dry and fire-prone 721 environments (Table 1). Given those dynamic changes in niche occupation, FE and FDv are 722 high in mountainous areas and in the semi-arid southern Mediterranean areas. Where T-NL 723 dominate FDv is slightly lower, i.e. niche occupation being less diverse. 724 5. Conclusion 725 We introduce a new generation, large-scale vegetation model for Europe allowing to 726 dynamically simulate functional diversity on a continental scale, alongside with productivity 727 and tree demography. Approximating competition of individual trees with randomly selected 728 functional trait combinations has proven successful to reproduce trait patterns, productivity 729 and tree demography of natural forest ecosystems. Whereas some trait ranges are 730 constrained by PFTs (SLA), others follow similar climatic gradients across all PFTs (WD). 731 These results complement the current spatially non-contiguous data of trait distributions 732 derived from earth observation data, and allow large-scale estimates of functional diversity 733 across Pan-Europe. We demonstrate how functional richness, divergence and evenness 734 vary strongly across Pan-Europe and emerge from environmental and competitive filtering 735 alike. Co-existence of functionally diverse trees results from plant trait diversity, tree 736 demography and disturbance under varying strength of environmental and competitive 737 filtering. 738 6. ACKNOWLEDGMENTS 739 This project has received funding from the European Union's Horizon 2020 research and 740 innovation programme under grant agreement No 641762 (ECOPOTENTIAL). The study 741 has been supported by the TRY initiative on plant traits (http://www.try-db.org). The TRY 742 initiative and database is hosted, developed and maintained by J. Kattge and G. Boenisch 743 (MaxPlanck Institute for Biogeochemistry, Jena, Germany) in collaboration with the German 744 Centre for Integrative Biodiversity Research (iDiv). We thank Giovanni Gligora for 745 contributing his data to TRY. TRY is been supported by DIVERSITAS core project 746 bioDISCOVERY, the IGBP, the Global Land Project, the UK Natural Environment Research 747 Council (NERC) through its program QUEST (Quantifying and Understanding the Earth 748 System) and the French programs 'Climat, Environnement, Societe' and 'Fondation pour la 749 Recherche sur la Biodiversite'. The contribution of M. Schaepman is supported by the 750 University of Zurich Research Priority Program on Global Change and Biodiversity (URPP 751 GCB). The research carried out at the Jet Propulsion Laboratory, California Institute of 19 752 Technology, was under a contract with the National Aeronautics and Space Administration. 753 Government sponsorship is acknowledged. 20 754 Figures 755 756 Figure 1 Conceptual scheme illustrating how flexible individual parameter are assigned to tree 757 individuals belonging to a specific PFT, including LES (leaf–) and SES (stem–economics) plant traits in 758 LPJmL-FIT. Parameter sets are assigned to individual trees at model initialization as well as gap opening. 759 Climate, soil conditions and competition between individual trees within and across PFTs (environmental 760 filtering) result in aggregated grid-cell trait distributions within patches. They are visualized as trait maps 761 and are used to quantify diversity indices. 762 763 764 Figure 2 Trait distribution of Specific Leaf Area (top panel) and wood density (bottom panel) for 765 broadleaved summergreen PFT (BL-S, left column), broadleaved evergreen PFT (BL-E, central column) 766 and both, boreal and temperate, needleleaved evergreen PFTs (NL, right column) as simulated by LPJmL- 767 FIT (light red) and observed by TRY data base entries (light cyan). Dashed vertical lines show the 768 simulated (light red) and measured (light cyan) average trait values for each PFT. 21 769 770 Figure 3 Mean Specific Leaf Area (SLA) and Wood Density (WD) distribution for each PFT plotted against 771 mean annual temperature (MAT) and mean annual precipitation (MAP). The colour of each dot represents 772 a mean SLA (a) and WD (b) value and each panel belongs to a PFT. (Broadleaved summergreen: BL-S; 773 Broadleaved evergreen: BL-E; Boreal needle-leaved evergreen: B-NL; and Temperate needle-leaved 774 evergreen: T-NL) Mean trait values are averaged over the 2004-2013 time period for 1000 patches per grid 775 cell. 776 777 Figure 4 Distribution of simulated trees in a 3-dimensional trait space (Height, WD, and SLA). We plot 778 density (far left) and functional diversity indices for 2’500 simulated patches (FR center left; FDv center 779 right; and FE (far right) for three sites (Seitseminen (top row), Lägern (middle row) and Dundo (bottom 780 row)). Ranges of the SLA and WD axes correspond to the maximum trait range across all PFTs used in 781 the simulation (see Table S2). Seitseminen, Laegern and Dundo represent boreal, mixed temperate and 782 Mediterranean-type forests, respectively. Each dot represents a trait combination of one tree larger than 783 5m while the colour indicates its PFT type. Plant Functional Types are: Broadleaved summergreen (BL- 784 S), Broadleaved Evergreen (BL-E), boreal needle-leaved evergreen (B-NL); and temperate needle-leaved 785 evergreen (T-NL). 22 786 787 Figure 5 Maps of functional richness (FR, panel a)), functional diversity (FDv, panel b)) and functional 788 evenness (FE, c)). All three functional diversity indices are based on SLA, LL, WD and tree height 789 averaged over all individual trees, which are at least 5m tall. Effects on FR, FDv and FE result from trait 790 variability within and across PFTs present in 1000 patches in a grid cell. 791 23 792 Tables PFT Resistance to Leaf strategy Growth strategy Environment (SLA range, Cold Warm Drought Fire SLA-LL temperature temperature relation) low medium low low Mostly thin and flexible WD, very flexible leaf strategy, but BL-S short-lived low WD in north, less specialized leaves high WD in south flexible growing strategy very competitive in temperate zone low high high medium moderate leaf high WD in dry less flexible leaf strategy BL-E cost environment conservative growing strategy intermediate lower WD under competitive in warm temperate leaf longevity wet conditions zones high low low low expensive and fast growth, due to conservative leaf strategy B-NL long-lived lower WD fast growing strategy needles specialist in cold climates with short growing season low high high high expensive and high WD in dry conservative leaf strategy T-NL long-lived environment, conservative growing strategy needles fast-growing under specialist in warm, dry and fire- colder climate, and prone environments thus lower WD 793 Table 1 Overview of the PFT characteristics: phenology, fire sensitivity, leaf and growth strategy and environment 1 794 795 7. REFERENCES 796 Bernard-Verdier, M., Navas, M.-L., Vellend, M., Violle, C., Fayolle, A., Garnier, E. & 797 Cornelissen, H. (2012) Community assembly along a soil depth gradient: contrasting 798 patterns of plant trait convergence and divergence in a Mediterranean rangeland. 799 Journal of Ecology, 100, 1422-1433. 800 Bohn, U., Zazanashvili, N. & Nakhutsrishvili, G. (2007) The Map of the Natural Vegetation of 801 Europe and its application in the Caucasus Ecoregion. Bulletin of the Georgia 802 Academy of Science, 175, 112-121. 803 Bruelheide, H., Dengler, J., Purschke, O., Lenoir, J., Jimenez-Alfaro, B., Hennekens, S.M., 804 Botta-Dukat, Z., Chytry, M., Field, R., Jansen, F., Kattge, J., Pillar, V.D., Schrodt, F., 805 Mahecha, M.D., Peet, R.K., Sandel, B., van Bodegom, P., Altman, J., Alvarez-Davila, 806 E., Arfin Khan, M.A.S., Attorre, F., Aubin, I., Baraloto, C., Barroso, J.G., Bauters, M., 807 Bergmeier, E., Biurrun, I., Bjorkman, A.D., Blonder, B., Carni, A., Cayuela, L., Cerny, 808 T., Cornelissen, J.H.C., Craven, D., Dainese, M., Derroire, G., De Sanctis, M., Diaz, 809 S., Dolezal, J., Farfan-Rios, W., Feldpausch, T.R., Fenton, N.J., Garnier, E., Guerin, 810 G.R., Gutierrez, A.G., Haider, S., Hattab, T., Henry, G., Herault, B., Higuchi, P., 811 Holzel, N., Homeier, J., Jentsch, A., Jurgens, N., Kacki, Z., Karger, D.N., Kessler, M., 812 Kleyer, M., Knollova, I., Korolyuk, A.Y., Kuhn, I., Laughlin, D.C., Lens, F., Loos, J., 813 Louault, F., Lyubenova, M.I., Malhi, Y., Marceno, C., Mencuccini, M., Muller, J.V., 814 Munzinger, J., Myers-Smith, I.H., Neill, D.A., Niinemets, U., Orwin, K.H., Ozinga, 815 W.A., Penuelas, J., Perez-Haase, A., Petrik, P., Phillips, O.L., Partel, M., Reich, P.B., 816 Romermann, C., Rodrigues, A.V., Sabatini, F.M., Sardans, J., Schmidt, M., Seidler, 817 G., Silva Espejo, J.E., Silveira, M., Smyth, A., Sporbert, M., Svenning, J.C., Tang, Z., 818 Thomas, R., Tsiripidis, I., Vassilev, K., Violle, C., Virtanen, R., Weiher, E., Welk, E., 819 Wesche, K., Winter, M., Wirth, C. & Jandt, U. (2018) Global trait-environment 820 relationships of plant communities. Nat Ecol Evol, 2, 1906-1917. 821 Butler, E.E., Datta, A., Flores-Moreno, H., Chen, M., Wythers, K.R., Fazayeli, F., Banerjee, 822 A., Atkin, O.K., Kattge, J., Amiaud, B., Blonder, B., Boenisch, G., Bond-Lamberty, B., 823 Brown, K.A., Byun, C., Campetella, G., Cerabolini, B.E.L., Cornelissen, J.H.C., 824 Craine, J.M., Craven, D., de Vries, F.T., Diaz, S., Domingues, T.F., Forey, E., 825 Gonzalez-Melo, A., Gross, N., Han, W., Hattingh, W.N., Hickler, T., Jansen, S., 826 Kramer, K., Kraft, N.J.B., Kurokawa, H., Laughlin, D.C., Meir, P., Minden, V., 827 Niinemets, U., Onoda, Y., Penuelas, J., Read, Q., Sack, L., Schamp, B., 828 Soudzilovskaia, N.A., Spasojevic, M.J., Sosinski, E., Thornton, P.E., Valladares, F., 829 van Bodegom, P.M., Williams, M., Wirth, C. & Reich, P.B. (2017) Mapping local and 830 global variability in plant trait distributions. Proc Natl Acad Sci U S A, 114, E10937- 831 E10946. 832 Cadotte, M.W., Carscadden, K. & Mirotchnick, N. (2011) Beyond species: functional diversity 833 and the maintenance of ecological processes and services. Journal of Applied 834 Ecology, 48, 1079-1087. 835 Carmona, C.P., de Bello, F., Mason, N.W.H. & Leps, J. (2016) Traits Without Borders: 836 Integrating Functional Diversity Across Scales. Trends Ecol Evol, 31, 382-394. 837 Carnicer, J., Barbeta, A., Sperlich, D., Coll, M. & Penuelas, J. (2013) Contrasting trait 838 syndromes in angiosperms and conifers are associated with different responses of 839 tree growth to temperature on a large scale. Front Plant Sci, 4, 409. 1 840 Chave, J., Coomes, D., Jansen, S., Lewis, S.L., Swenson, N.G. & Zanne, A.E. (2009) 841 Towards a worldwide wood economics spectrum. Ecol Lett, 12, 351-66. 842 D'Odorico, P., Gonsamo, A., Pinty, B., Gobron, N., Coops, N., Mendez, E., & Schaepman, 843 M.E. (2014) Intercomparison of fraction of absorbed photosynthetically active 844 radiation products derived from satellite data over Europe. Remote Sensing of 845 Environment, 142, 141-154. 846 Díaz, S., Kattge, J., Cornelissen, J.H.C., Wright, I.J., Lavorel, S., Dray, S., Reu, B., Kleyer, 847 M., Wirth, C., Colin Prentice, I., Garnier, E., Bönisch, G., Westoby, M., Poorter, H., 848 Reich, P.B., Moles, A.T., Dickie, J., Gillison, A.N., Zanne, A.E., Chave, J., Joseph 849 Wright, S., Sheremet’ev, S.N., Jactel, H., Baraloto, C., Cerabolini, B., Pierce, S., 850 Shipley, B., Kirkup, D., Casanoves, F., Joswig, J.S., Günther, A., Falczuk, V., Rüger, 851 N., Mahecha, M.D. & Gorné, L.D. (2016) The global spectrum of plant form and 852 function. Nature, 529, 167-171. 853 Douma, J.C., Witte, J.-P.M., Aerts, R., Bartholomeus, R.P., Ordoñez, J.C., Venterink, H.O., 854 Wassen, M.J. & van Bodegom, P.M. (2012) Towards a functional basis for predicting 855 vegetation patterns; incorporating plant traits in habitat distribution models. 856 Ecography, 35, 294-305. 857 Ellis, E.C., Kaplan, J.O., Fuller, D.Q., Vavrus, S., Goldewijk, K.K. & Verburg, P.H. (2013) 858 Used planet: A global history. Proceedings of the National Academy of Sciences of 859 the United States of America, 110, 7978-7985. 860 FAO (2018) Terms and Definitions FRA 2020. In: Forest Resources Assessment Working 861 Papers (ed. D.P.B.A. Pekkarinen.). FAO Forestry Department, Rome. 862 Forkel, M., Carvalhais, N., Schaphoff, S., von Bloh, W., Migliavacca, M., Thurner, M. & 863 Thonicke, K. (2014) Identifying environmental controls on vegetation greenness 864 phenology through model-data integration. Biogeosciences, 11, 7025-7050. 865 Funk, J.L., Larson, J.E., Ames, G.M., Butterfield, B.J., Cavender-Bares, J., Firn, J., Laughlin, 866 D.C., Sutton-Grier, A.E., Williams, L. & Wright, J. (2017) Revisiting the Holy Grail: 867 using plant functional traits to understand ecological processes. Biol Rev Camb 868 Philos Soc, 92, 1156-1173. 869 Garnier, E., Navas, M.-L. & Grigulis, K. (2016) Plant Functional Diversity. Organisms, traits, 870 community structure, and ecosystem properties. Oxford University Press, Oxford. 871 Garonna, I., de Jong, R., Stöckli, R., Schmid, B., Schenkel, D., David, S., & Schaepman, 872 M.E. (2018) Shifting relative importance of climatic constraints on land surface 873 phenology. Environmental Research Letters, 13, 024025. 874 Geller, G.N., Halpin, P.N., Helmuth, B., Hestir, E.L., Skidmore, A., Abrams, M.J., Aguirre, N., 875 Blair, M., Botha, E., Collof, M., Dawson, T., Franklin, J., Horning, N., James, C., 876 Magnusson, W., Santos, M.J., Schill, S.R. & Williams, K. (2017) Remote Sensing for 877 Biodiversity. The GEO Handbook on Biodiversity Observation Networks (ed. by M. 878 Walters and R.J. Scholes), pp. 187-210. Springer Nature, Cham, Switzerland. 879 Hallik, L., Niinemets, Ü. & Wright, I. (2009) Are species shade and drought tolerance 880 reflected in leaf-level structural and functional differentiation in Northern Hemisphere 881 temperate woody flora? . The New Phytologist, 184, 257-274. 882 Healey, S.P., Hernandez, M.W., Edwards, D.P., Lefsky, M.A., Freeman, E., Patterson, P.L., 883 Lindquist, E.J. & Lister, A.J. (2015) CMS: GLAS LiDAR-derived Global Estimates of 884 Forest Canopy Height, 2004-2008. In: (ed. O. Daac), Oak Ridge, Tennessee, USA. 885 Hisano, M., Searle, E.B. & Chen, H.Y.H. (2018) Biodiversity as a solution to mitigate climate 886 change impacts on the functioning of forest ecosystems. Biological Reviews, 93, 439- 887 456. 2 888 Hooper, D.U., Chapin, F.S., Ewel, J.J., Hector, A., Inchausti, P., Lavorel, S., Lawton, J.H., 889 Lodge, D.M., Loreau, M., Naeem, S., Schmid, B., Setala, H., Symstad, A.J., 890 Vandermeer, J. & Wardle, D.A. (2005) EFFECTS OF BIODIVERSITY ON 891 ECOSYSTEM FUNCTIONING: A CONSENSUS OF CURRENT KNOWLEDGE. 892 Ecological Monographs, 75, 3-35. 893 Jetz, W., Cavender-Bares, J., Pavlick, R., Schimel, D., Davis, F.W., Asner, G.P., Guralnick, 894 R., Kattge, J., Latimer, A.M., Moorcroft, P., Schaepman, M.E., Schildhauer, M.P., 895 Schneider, F.D., Schrodt, F., Stahl, U. & Ustin, S.L. (2016) Monitoring plant functional 896 diversity from space. Nat Plants, 2, 16024. 897 Kattge, J., S. Díaz, S. Lavorel, I. C. Prentice, P. Leadley & G. Bönisch, E.G., M. Westoby, P. 898 B. Reich, I. J. Wright, J. H. C. Cornelissen, C. Violle, Harrison, S. P., , P. M. Van 899 Bodegom, M. Reichstein, B. J. Enquist, N. A. Soudzilovskaia, D. D. Ackerly, M. 900 Anand, O. Atkin, M. Bahn, T. R. Baker, D. Baldocchi, R. Bekker, C. C. Blanco, B. 901 Blonder, W. J. Bond, R. Bradstock, D. E. Bunker, F. Casanoves, J. Cavender-Bares, 902 J. Q. Chambers, F. S. Chapin Iii, J. Chave, D. Coomes, W. K. Cornwell, J. M. Craine, 903 B. H. Dobrin, L. Duarte, W. Durka, J. Elser, G. Esser, M. Estiarte, W. F. Fagan, J. 904 Fang, F. Fernández-Méndez, A. Fidelis, And B. Finegan, Flores, O., Ford, H., Frank, 905 D., Freschet, G. T., Fyllas, N. M., Gallagher, R. V., Green, W. A., Gutierrez, A. G., 906 Hickler, T., Higgins, S. I., Hodgson, J. G., Jalili, A., Jansen, S., Joly, C. A., Kerkhoff, 907 A. J., Kirkup, D., Kitajima, K., Kleyer, M., Klotz, S., Knops, J. M. H., Kramer, K., 908 Kühn, I., Kurokawa, H., Laughlin, D., Lee, T. D., Leishman, M., Lens, F., Lenz, T., 909 Lewis, S. L., Lloyd, J., Llusià, J., Louault, F., Ma, S., Mahecha, M. D., Manning, P., 910 Massad, T., Medlyn, B. E., Messier, J., Moles, A. T., Müller, S. C., Nadrowski, K., 911 Naeem, S., Niinemets, Ü., Nöllert, S., Nüske, A., Ogaya, R., Oleksyn, J., 912 Onipchenko, V. G., Onoda, Y., Ordoñez, J., Overbeck, G., Ozinga, W. A., Patiño, S., 913 Paula, S., Pausas, J. G., Peñuelas, J., Phillips, O. L., Pillar, V., Poorter, H., Poorter, 914 L., Poschlod, P., Prinzing, A., Proulx, R., Rammig, A., Reinsch, S., Reu, B., Sack, L., 915 Salgado-Negret, B., Sardans, J., Shiodera, S., Shipley, B., Siefert, A., Sosinski, E., 916 Soussana, J.-F., Swaine, E., Swenson, N., Thompson, K., Thornton, P., Waldram, 917 M., Weiher, E., White, M., White, S., Wright, S. J., Yguel, B., Zaehle, S., Zanne, A. E. 918 and Wirth, C (2011) TRY - a global database of plant traits. Global Change Biology, 919 17, 2905–2935. 920 King, D.A., Davies, S. J., Tan, S., & Noor, N. S. M. (2006) The role of wood density and stem 921 support costs in the growth and mortality of tropical trees. Journal of Ecology, 94, 922 670-680. 923 Kunstler, G., Falster, D., Coomes, D.A., Hui, F., Kooyman, R.M., Laughlin, D.C., Poorter, L., 924 Vanderwel, M., Vieilledent, G., Wright, S.J., Aiba, M., Baraloto, C., Caspersen, J., 925 Cornelissen, J.H., Gourlet-Fleury, S., Hanewinkel, M., Herault, B., Kattge, J., 926 Kurokawa, H., Onoda, Y., Penuelas, J., Poorter, H., Uriarte, M., Richardson, S., Ruiz- 927 Benito, P., Sun, I.F., Stahl, G., Swenson, N.G., Thompson, J., Westerlund, B., Wirth, 928 C., Zavala, M.A., Zeng, H., Zimmerman, J.K., Zimmermann, N.E. & Westoby, M. 929 (2016) Plant functional traits have globally consistent effects on competition. Nature, 930 529, 204-7. 931 Laanisto, L. & Niinemets, Ü. (2015) Polytolerance to abiotic stresses: how universal is the 932 shade-drought tolerance trade-off in woody species? Global Ecology & 933 Biogeography, 24, 571-580. 934 Laliberté, E., Legendre, P., & Shipley, B. (2014) FD: Measuring functional diversity (FD) from 935 multiple traits, and other tools for functional ecology. R package version 1.0-12. 936 Langan, L., Higgins, S.I. & Scheiter, S. (2017) Climate‐biomes, pedo‐biomes or pyro‐biomes: 937 which world view explains the tropical forest–savanna boundary in South America? 938 Journal of Biogeography, 44, 2319-2330. 3 939 Lausch, A., Bastian, O., Klotz, S., Leitão, P.J., Jung, A., Rocchini, D., Schaepman, M.E., 940 Skidmore, A.K., Tischendorf, L., & Knapp, S. (2018) Understanding and assessing 941 vegetation health by in situ species and remote-sensing approaches. Methods in 942 Ecology and Evolution, 9, 1799-1809. 943 Lavorel, S. & Garnier, E. (2002a) Predicting Changes in Community Composition and 944 Ecosystem Functioning from Plant Traits: 945 Revisiting the Holy Grail. Functional Ecology, 16, 545-556. 946 Lavorel, S. & Garnier, E. (2002b) Predicting changes in community composition and 947 ecosystem functioning from plant traits: revisiting the Holy Grail. Functional Ecology, 948 16, 545-556. 949 Lavorel, S., Diaz, S., Cornelissen, J.H.C., Garnier, E., Harrison, S.P., McIntyre, S., Pausas, 950 J.G., Perez-Harguinedeguy, Roumet, C. & Urcelay, C. (2007) Plant Functional 951 Types: Are We Getting Any Closer to the Holy Grail? Terrestrial Ecosystems in a 952 Changing World (ed. by C. Jp, P. D and P. Lf), pp. 149-160. Springer, Berlin, 953 Heidelberg, New York. 954 Lefsky, M.A., M. Keller, Y. Pang, P.B. De Camargo, and M.O. Hunter (2007) Revised 955 method for forest canopy height estimation from Geoscience Laser Altimeter System 956 waveforms. Journal of Applied Remote Sensing, 1, 013537-013518. 957 Liang, J., Crowther, T.W., Picard, N., Wiser, S., Zhou, M., Alberti, G., Schulze, E.D., 958 McGuire, A.D., Bozzato, F., Pretzsch, H., de-Miguel, S., Paquette, A., Herault, B., 959 Scherer-Lorenzen, M., Barrett, C.B., Glick, H.B., Hengeveld, G.M., Nabuurs, G.J., 960 Pfautsch, S., Viana, H., Vibrans, A.C., Ammer, C., Schall, P., Verbyla, D., 961 Tchebakova, N., Fischer, M., Watson, J.V., Chen, H.Y., Lei, X., Schelhaas, M.J., Lu, 962 H., Gianelle, D., Parfenova, E.I., Salas, C., Lee, E., Lee, B., Kim, H.S., Bruelheide, 963 H., Coomes, D.A., Piotto, D., Sunderland, T., Schmid, B., Gourlet-Fleury, S., Sonke, 964 B., Tavani, R., Zhu, J., Brandl, S., Vayreda, J., Kitahara, F., Searle, E.B., Neldner, 965 V.J., Ngugi, M.R., Baraloto, C., Frizzera, L., Balazy, R., Oleksyn, J., Zawila- 966 Niedzwiecki, T., Bouriaud, O., Bussotti, F., Finer, L., Jaroszewicz, B., Jucker, T., 967 Valladares, F., Jagodzinski, A.M., Peri, P.L., Gonmadje, C., Marthy, W., O'Brien, T., 968 Martin, E.H., Marshall, A.R., Rovero, F., Bitariho, R., Niklaus, P.A., Alvarez-Loayza, 969 P., Chamuya, N., Valencia, R., Mortier, F., Wortel, V., Engone-Obiang, N.L., Ferreira, 970 L.V., Odeke, D.E., Vasquez, R.M., Lewis, S.L. & Reich, P.B. (2016) Positive 971 biodiversity-productivity relationship predominant in global forests. Science, 354 972 Liebergesell, M., Reu, B., Stahl, U., Freiberg, M., Welk, E., Kattge, J., Cornelissen, J.H., 973 Penuelas, J. & Wirth, C. (2016) Functional Resilience against Climate-Driven 974 Extinctions - Comparing the Functional Diversity of European and North American 975 Tree Floras. PLoS One, 11, e0148607. 976 Madrigal-Gonzalez, J., Ruiz-Benito, P., Ratcliffe, S., Calatayud, J., Kandler, G., Lehtonen, 977 A., Dahlgren, J., Wirth, C. & Zavala, M.A. (2016) Complementarity effects on tree 978 growth are contingent on tree size and climatic conditions across Europe. Scientific 979 Reports, 6, 10. 980 Mason, N.W.H. & de Bello, F. (2013) Functional diversity: a tool for answering challenging 981 ecological questions. Journal of Vegetation Science, 24, 777-780. 982 Mason, N.W.H., Mouillot, D., Lee, W.G. & Wilson, J.B. (2005) Functional richness, functional 983 evenness and functional divergence: the primary components of functional diversity. 984 Oikos, 111, 112-118. 985 Mason, N.W.H., de Bello, F., Mouillot, D., Pavoine, S. & Dray, S. (2013) A guide for using 986 functional diversity indices to reveal changes in assembly processes along ecological 987 gradients. Journal of Vegetation Science, 24, 794-806. 4 988 Mitchell, T.D. (2002) A comprehensive set of climate scenarios for Europe. in prep., 989 Mori, A.S., Lertzman, K.P. & Gustafsson, L. (2017) Biodiversity and ecosystem services in 990 forest ecosystems: a research agenda for applied forest ecology. Journal of Applied 991 Ecology, 54, 12-27. 992 Mouillot, D., Graham, N.A.J., Villeger, S., Mason, N.W.H. & Bellwood, D.R. (2013) A 993 functional approach reveals community responses to disturbances. Trends in 994 Ecology & Evolution, 28, 167-177. 995 Nachtergaele, F., van Velthuizen, H., Verelst, L., Batjes, N., Dijkshoorn, K., van Engelen, V., 996 Fischer, G., Jones, A., & Montanarella, L., and Petri, M. (2009) Harmonized world 997 soil database, Food and Agriculture Organization of the United Nations. In: (ed. Fao), 998 Rome. 999 Naeem, S., Duffy, J.E. & Zavaleta, E.S. (2012) The functions of biological diversity in an Age 1000 of Extinction. Science, 336, 1401-1406. 1001 Naveh, Z. (1990) Fire in the Mediterranean - a landscape ecological perspective. Fire in 1002 ecosystem dynamics: Mediterranean and Northern perspectives (ed. by J.G. 1003 Goldammer and M.J. Jenkins). SPB Academic Publ., The Hague. 1004 Niinemets, Ü. & Valladares, F. (2006) Tolerance to shade, drought and waterlogging in the 1005 temperate dendroflora of the Northern hemisphere: tradeoffs, phylogenetic signal and 1006 implications for niche differentiation. Ecological Monographs, 76, 521-547. 1007 Pakeman, R.J. (2011) Functional diversity indices reveal the impacts of land use 1008 intensification on plantcommunity assembly. Journal of Ecology, 99, 1143-1151. 1009 Papale, D., Reichstein, M., Aubinet, M., Canfora, E., Bernhofer, C., Kutsch, W. & Yakir, D. 1010 (2006) Towards a standardized processing of Net Ecosystem Exchange measured 1011 with eddy covariance technique: algorithms and uncertainty estimation. 1012 Biogeosciences, 3, 571-583. 1013 Paul-Limoges, E., Damm, A., Hueni, A., Liebisch, F., Eugster, W., Schaepman, M.E., & 1014 Buchmann, N. (2018) Effect of environmental conditions on sun-induced 1015 fluorescence in a mixed forest and a cropland. Remote Sensing of Environment, 219, 1016 310-323. 1017 Prentice, I., Bondeau, A., Cramer, W., Harrison, S., Hickler, T., Lucht, W., Sitch, S., Smith, 1018 B. & Sykes, M. (2007) Dynamic global vegetation modelling: quantifying terrestrial 1019 ecosystem responses to large-scale environmental change. . Terrestrial Ecosystems 1020 in a Changing World (ed. by C. Jp, P. D and P. Lf). Springer, Berlin, Heidelberg, New 1021 York. 1022 Ratcliffe, S., Liebergesell, M., Ruiz-Benito, P., Madrigal González, J., Muñoz Castañeda, 1023 J.M., Kändler, G., Lehtonen, A., Dahlgren, J., Kattge, J., Peñuelas, J., Zavala, M.A. & 1024 Wirth, C. (2016) Modes of functional biodiversity control on tree productivity across 1025 the European continent. Global Ecology and Biogeography, 25, 251-262. 1026 Ratcliffe, S., Wirth, C., Jucker, T., van der Plas, F., Scherer-Lorenzen, M., Verheyen, K., 1027 Allan, E., Benavides, R., Bruelheide, H., Ohse, B., Paquette, A., Ampoorter, E., 1028 Bastias, C.C., Bauhus, J., Bonal, D., Bouriaud, O., Bussotti, F., Carnol, M., 1029 Castagneyrol, B., Checko, E., Dawud, S.M., De Wandeler, H., Domisch, T., Finer, L., 1030 Fischer, M., Fotelli, M., Gessler, A., Granier, A., Grossiord, C., Guyot, V., Haase, J., 1031 Hattenschwiler, S., Jactel, H., Jaroszewicz, B., Joly, F.X., Kambach, S., Kolb, S., 1032 Koricheva, J., Liebersgesell, M., Milligan, H., Muller, S., Muys, B., Nguyen, D., Nock, 1033 C., Pollastrini, M., Purschke, O., Radoglou, K., Raulund-Rasmussen, K., Roger, F., 1034 Ruiz-Benito, P., Seidl, R., Selvi, F., Seiferling, I., Stenlid, J., Valladares, F., Vesterdal, 1035 L. & Baeten, L. (2017) Biodiversity and ecosystem functioning relations in European 1036 forests depend on environmental context. Ecology Letters, 20, 1414-1426. 5 1037 Reich, P.B. (2014) The world-wide ‘fast–slow’ plant economics spectrum: a traits manifesto. 1038 Journal of Ecology, 102, 275-301. 1039 Reichstein, M., Falge, E., Baldocchi, D., Papale, D., Aubinet, M., Berbigier, P. & Grünwald, 1040 T. (2005) On the separation of net ecosystem exchange into assimilation and 1041 ecosystem respiration: review and improved algorithm. Global Change Biology, 11, 1042 1424-1439. 1043 Ruiz-Benito, P., Ratcliffe, S., Jump, A.S., Gomez-Aparicio, L., Madrigal-Gonzalez, J., Wirth, 1044 C., Kandler, G., Lehtonen, A., Dahlgren, J., Kattge, J. & Zavala, M.A. (2017) 1045 Functional diversity underlies demographic responses to environmental variation in 1046 European forests. Global Ecology and Biogeography, 26, 128-141. 1047 Running, S. (2015) MODIS/TERRA Gross Primary Productivity 8-Day L4 Global 500 m SIN 1048 Grid V006. In: (ed. U.E.R.O.a.S.E.C. Nasa Eosdis Land Processes Daac), Sioux 1049 Falls, South Dakota (Https://Lpdaac.Usgs.Gov). 1050 Sakschewski, B., von Bloh, W., Boit, A., Rammig, A., Kattge, J., Poorter, L., Penuelas, J. & 1051 Thonicke, K. (2015) Leaf and stem economics spectra drive diversity of functional 1052 plant traits in a dynamic global vegetation model. Global Change Biology, 21, 2711- 1053 2725. 1054 Sandel, B., Gutierrez, A.G., Reich, P.B., Schrodt, F., Dickie, J. & Kattge, J. (2015) Estimating 1055 the missing species bias in plant trait measurements. Journal of Vegetation Science, 1056 26, 828-838. 1057 Schaphoff, S., Forkel, M., Müller, C., Knauer, J., von Bloh, W., Gerten, D., Jägermeyr, J., 1058 Lucht, W., Rammig, A., Thonicke, K. & Waha, K. (2018a) LPJmL4 – a dynamic global 1059 vegetation model with managed land – Part 2: Model evaluation. Geoscientific Model 1060 Development, 11, 1377-1403. 1061 Schaphoff, S., von Bloh, W., Rammig, A., Thonicke, K., Biemans, H., Forkel, M., Gerten, D., 1062 Heinke, J., Jägermeyr, J., Knauer, J., Langerwisch, F., Lucht, W., Müller, C., Rolinski, 1063 S. & Waha, K. (2018b) LPJmL4 – a dynamic global vegetation model with managed 1064 land – Part 1: Model description. Geoscientific Model Development, 11, 1343-1375. 1065 Schneider, F.D., Morsdorf, F., Schmid, B., Petchey, O.L., Hueni, A., Schimel, D.S. & 1066 Schaepman, M.E. (2017) Mapping functional diversity from remotely sensed 1067 morphological and physiological forest traits. Nat Commun, 8, 1441. 1068 Schneider, U., Becker, A., Finger, P., Meyer-Christoffer, A., Rudolf, B. & Ziese, M. (2011) 1069 GPCC Full Data Reanalysis Version 6.0 at 0.5°: Monthly Land-Surface Precipitation 1070 from Rain-Gauges built on GTS-based and Historic Data. 1071 Šímová, I., Violle, C., Svenning, J.-C., Kattge, J., Engemann, K., Sandel, B., Peet, R.K., 1072 Wiser, S.K., Blonder, B., McGill, B.J., Boyle, B., Morueta-Holme, N., Kraft, N.J.B., 1073 van Bodegom, P.M., Gutiérrez, A.G., Bahn, M., Ozinga, W.A., Tószögyová, A. & 1074 Enquist, B.J. (2018) Spatial patterns and climate relationships of major plant traits in 1075 the New World differ between woody and herbaceous species. Journal of 1076 Biogeography, 45, 895-916. 1077 Sitch, S., Smith, B., Prentice, I.C., Arneth, A., Bondeau, A., Cramer, W., Kaplan, J.O., Levis, 1078 S., Lucht, W., Sykes, M.T., Thonicke, K. & Venevsky, S. (2003) Evaluation of 1079 ecosystem dynamics, plant geography and terrestrial carbon cycling in the LPJ 1080 Dynamic Global Vegetation Model. Global Change Biology, 9, 161-185. 1081 Sitch, S., Huntingford, C., Gedney, N., Levy, P.E., Lomas, M., Piao, S.L., Betts, R.A., Ciais, 1082 P., Cox, P.M., Friedlingstein, P., Jones, C.D., Prentice, I.C. & Woodward, F.I. (2008) 1083 Evaluation of the terrestrial carbon cycle, future plant geography and climate-carbon 1084 cycle feedbacks using five Dynamic Global Vegetation Models (DGVMs). Global 1085 Change Biology, 14, 2015-2039. 6 1086 Thonicke, K., Venevsky, S., Sitch, S. & Cramer, W. (2001) The role of fire disturbance for 1087 global vegetation dynamics: coupling fire into a Dynamic Global Vegetation Model. 1088 Global Ecology and Biogeography, 10, 661-677. 1089 Thurner, M., Beer, C., Santoro, M., Carvalhais, N., Wutzler, T., Schepaschenko, D., ... & 1090 Schmullius, C. (2014) Carbon stock and density of northern boreal and temperate 1091 forests. Global Ecology and Biogeography, 23, 297-310. 1092 Tinner, W., Hubschmid, P., Wehrli, M., Ammann, B., & Conedera, M. (1999) Long‐term 1093 forest fire ecology and dynamics in southern Switzerland. Journal of Ecology, 87, 1094 273-289. 1095 Van Bodegom, P.M., Douma, J.C., Witte, J.P.M., Ordoñez, J.C., Bartholomeus, R.P. & 1096 Aerts, R. (2012) Going beyond limitations of plant functional types when predicting 1097 global ecosystem-atmosphere fluxes: exploring the merits of traits-based 1098 approaches. Global Ecology and Biogeography, 21, 625-636. 1099 Villeger, S., Mason, N.W.H. & Mouillot, D. (2008) New multidimensional functional diversity 1100 indices for a multifaceted framework in functional ecology. Ecology, 89, 2290-2301. 1101 Weedon, G.P., Balsamo, G., Bellouin, N., Gomes, S., Best, M. J., & Viterbo, P. (2014) The 1102 WFDEI meteorological forcing data set: WATCH Forcing Data methodology applied 1103 to ERA-Interim reanalysis data Water Resources Research, 50, 7505–7514. 1104 Weedon, G.P., Gomes, S., Viterbo, P., Shuttleworth, W. J., Blyth, E., Österle, H., … Best, M. 1105 (2011) Creation of the WATCH Forcing Data and Its Use to Assess Global and 1106 Regional Reference Crop Evaporation over Land during the Twentieth Century. 1107 Journal of Hydrometeorology, 12 823–848. 1108 Wright, I.J., Ning Dong, Maire, V., Prentice, I.C., Westoby, M., Díaz, S., Gallagher, R.V., 1109 Jacobs, B.F., Kooyman, R., Law, E.A., Leishman, M.R., Niinemets, Ü., Reich, P.B., 1110 Sack, L., Villar, R., Wang, H. & Wilf, P. (2017) Global climatic drivers of leaf size. 1111 Science, 357, 917-921. 1112 Wright, I.J., Peter B. Reich, Mark Westoby, David D. Ackerly, Zdravko Baruch, Frans 1113 Bongers, Jeannine Cavender-Bares, Terry Chapin, Johannes H. C. Cornelissen, 1114 Matthias Diemer, Jaume Flexas, Eric Garnier, Philip K. Groom, Javier Gulias, Kouki 1115 Hikosaka, Byron B. Lamont, Tali Lee, William Lee, Christopher Lusk, Jeremy J. 1116 Midgley, Marie-Laure Navas, Ülo Niinemets, Jacek Oleksyn, Noriyuki Osada, Hendrik 1117 Poorter, Pieter Poot, Lynda Prior, Vladimir I. Pyankov, Catherine Roumet, Sean C. 1118 Thomas, Mark G. Tjoelker, Veneklaas, E.J. & Villar, R. (2004) The worldwide leaf 1119 economics spectrum. Nature 428, 821-827. 1120 Zhang, T., Niinemets, U., Sheffield, J. & Lichstein, J.W. (2018) Shifts in tree functional 1121 composition amplify the response of forest biomass to climate. Nature, 556, 99-102. 1122 7