Received: 6 May 2019 | Revised: 5 December 2019 | Accepted: 12 December 2019 DOI: 10.1111/jbi.13809 RE SE ARCH PAPER Simulating functional diversity of European natural forests along climatic gradients Kirsten Thonicke1 | Maik Billing1,2 | Werner von Bloh1 | Boris Sakschewski1 | Ülo Niinemets3 | Josep Peñuelas4,5 | J. Hans C. Cornelissen6 | Yusuke Onoda7 | Peter van Bodegom8 | Michael E. Schaepman9 | Fabian D. Schneider10 | Ariane Walz2 1 Research Domain 1 “Earth System Analysis”, Potsdam Institute for Climate Abstract Impact Research (PIK), Member of the Aim: We analyse how functional diversity (FD) varies across European natural forests Leibniz Association, Potsdam, Germany 2 to understand the effects of environmental and competitive filtering on plant trait Institute of Environmental Science and Geography, University of Potsdam, distribution. Potsdam-Golm, Germany Location: Forest ecosystems in Europe from 11°W to 36°E and 29.5°N to 62°N. 3 Estonian University of Life Sciences, Tartu, Estonia Taxon: Pinaceae, Fagaceae and Betulaceae, Oleaceae, Tiliaceae, Aceraceae, 4 National Research Council (CSIC), Global Leguminosae (unspecific). Ecology Unit at Center for Ecological Methods: We adopted the existing Dynamic Global Vegetation Model Lund-Potsdam- Research and Forestry Applications (CREAF), Universitat Autonoma de Jena managed Land of flexible individual traits (LPJmL-FIT) for Europe by eliminating Barcelona (UAB), Bellaterra, Catalonia, Spain both bioclimatic limits of plant functional types (PFTs) and replacing prescribed val- 5 Center for Ecological Research and ues of functional traits for PFTs with emergent values under influence of environ- Forestry Applications (CREAF), Bellaterra, Catalonia, Spain mental filtering and competition. We quantified functional richness (FR), functional 6 Systems Ecology, Department of Ecological divergence (FDv) and functional evenness (FE) in representative selected sites and at Science, Vrije Universiteit, Amsterdam, the Netherlands Pan-European scale resulting from simulated functional and structural trait combina- 7 Division of Environmental Science and tions of individual trees. While FR quantifies the amount of occupied trait space, FDv Technology, Graduate School of Agriculture, and FE describe the distribution and abundance of trait combinations, respectively, Kyoto University, Kyoto, Japan 8 Institute of Environmental Sciences, in a multidimensional trait space. Department Environmental Biology, Leiden Results: Lund-Potsdam-Jena managed Land of flexible individual traits reproduces University, Leiden, the Netherlands 9 spatial PFTs and local trait distributions and agrees well with observed productivity, Remote Sensing Laboratories, Dept. of Geography, University of Zurich, Zurich, biomass and tree height of European natural forests. The observed site-specific trait Switzerland distributions and spatial gradients of traits of the leaf- and stem-resource economics 10 Jet Propulsion Laboratory, California spectra coincide with environmental filtering and the competition for light and water Institute of Technology, Pasadena, CA, USA in environments with strong abiotic stress. Where deciduous and needle-leaved trees Correspondence co-occur, for example, in boreal and mountainous forests, the potential niche space is Kirsten Thonicke, Research Domain 1 “Earth System Analysis”, Potsdam Institute for wide (high FR), and extreme ends in the niche space are occupied (high FDv). We find Climate Impact Research (PIK), Member of high FDv in Mediterranean forests where drought increasingly limits tree growth, the Leibniz Association, Potsdam, Germany. Email: kirsten.thonicke@pik-potsdam.de thus niche differentiation becomes more important. FDv decreases in temperate for- ests where a cold climate increasingly limits growth efficiency of broad-leaved sum- mer green trees, thus reducing the importance of competitive exclusion. Highest FE This is an open access article under the terms of the Creat​ive Commo​ns Attri​bution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2020 The Authors. Journal of Biogeography published by John Wiley & Sons Ltd Journal of Biogeography. 2020;47:1069–1085. wileyonlinelibrary.com/journal/jbi | 1069 1070 | THONICKE et al. Funding information European Union's Horizon 2020 Research was simulated in wet Atlantic and southern Europe which indicated relatively even and Innovation Programme, Grant/Award niche occupation and thus high resource-use efficiency. Number: 641762; University of Zurich Research Priority Program on Global Main Conclusions: We find FD resulting from both environmental and competitive Change and Biodiversity filtering. Pan-European FR, FDv and FE demonstrate the influence of climate gradi- Handling Editor: Holger Kreft ents and intra- and inter-PFT competition. The indices underline a generally high FD of natural forests in Europe. Co-existence of functionally diverse trees across PFTs emerges from alternative (life-history) strategies, disturbance and tree demography. KEYWORDS boreal forest, co-existence, competitive filtering, dynamic global vegetation model, environmental filtering, flexible individual traits, mediterranean forest, natural forest, temperate forest 1 | I NTRO D U C TI O N Wright et al., 2004). To quantify how environmental and competitive filtering influence niche complementarity, multidimensional indices of Functional diversity (FD) is a key control of the stability and adapt- FD are required to quantify occupation and overlap of niches (Mason ability of ecosystems under climate change (Yachi & Loreau, 1999). & de Bello, 2013). In addition, these indices should be scale indepen- Abiotic conditions (e.g. climate, soil) as well as biotic processes (e.g. dent and applicable at regional scale to investigate changes in FD competition) determine plant community assembly (Kunstler et al., along climatic gradients (Carmona, Bello, Mason, & Leps, 2016). 2016; Ratcliffe et al., 2017; Ruiz-Benito et al., 2017), and thus the FD Many different indices have been used to describe the size of trait of communities (Cadotte, Carscadden, & Mirotchnick, 2011; Lavorel spaces, the distribution and the clustering of their trait combinations & Garnier, 2002a; Naeem, Duffy, & Zavaleta, 2012). Ecosystem FD (e.g. Ratcliffe et al., 2016; Schneider et al., 2017). Villeger et al. (2008) has been measured based on the diversity of morphological and/or suggested to describe these FD aspects using three independent in- physiological plant traits (Villeger, Mason, & Mouillot, 2008), which dices: functional richness (FR), divergence (FDv) and evenness (FE). are linked to plant productivity, transpiration or nutrient cycling, and Whereas FR quantifies the amount of occupied trait space, FDv and thus to ecosystem functions (cf. Lavorel & Garnier, 2002a). For in- FE describe the distribution and abundance of trait combinations in a stance, specific leaf area (SLA) and leaf nitrogen content have been multidimensional trait space (Mason et al., 2005; Villeger et al., 2008). linked to plant productivity, and stem traits such as wood density Whereas FR describes the size of potentially available, functional (WD) to carbon storage and evapotranspiration (Funk et al., 2017) space, in which niches can be occupied by plants, FDv quantifies and, among others, determine species tolerance of environmental the distribution of trait values, thus the degree of niche differentia- stresses (Hallik, Niinemets, & Wright, 2009) and control competitive tion, likely the result from competitive exclusion (Garnier, Navas, & interactions between individual plants (Kunstler et al., 2016). Grigulis, 2016; Mason et al., 2005). FE describes the regularity of trait Climatic conditions are the most important drivers of commu- distribution and points to resource-use efficiency within the occupied nity assembly and generally constrain the relations between traits trait space. Lower FR could relate to lower capability of an ecosystem globally (Butler et al., 2017; Šímová et al., 2018). Gradients in trait to buffer environmental stress, whereas lower FE and FDv could indi- expression, associated with climatic conditions, have been found cate reduced ecosystem resilience (Mason et al., 2005). at global scales (Díaz et al., 2016; Wright et al., 2017) and at the By dividing FD into richness, divergence and evenness, the mech- European scale by extrapolating site-specific plant trait data (Butler anisms that link biodiversity to ecosystem functions and describe et al., 2017) and using forest inventories (Ruiz-Benito et al., 2017). community assembly can be described (Mason, Bello, Mouillot, Recent studies investigated a mix of globally important physiolog- Pavoine, & Dray, 2013; Mason et al., 2005). Computing scale-inde- ical traits (e.g. SLA, WD, seed mass and leaf nitrogen content) and pendent FD indices from trait distributions (Carmona et al., 2016) morphological traits (including maximum plant height and basal allows investigating FD from community to meta-community scale area) (Ratcliffe et al., 2016; Ruiz-Benito et al., 2017), or separate and extrapolating them to the regional or continental scale. Huge ef- the effects between those trait types (e.g. Madrigal-Gonzalez et al., forts are under way to explore links among plant traits, vegetation 2016; Schneider et al., 2017). They advance our understanding of composition and climate based on site data (Bruelheide et al., 2018; spatial pattern of plant traits at the landscape and continental scale. Díaz et al., 2016; Wright et al., 2017). They are complemented by Under given climate conditions, the distribution of multiple in- dynamic global vegetation models (DGVM) with flexible or adaptive dividual traits results from community assembly rules and opens a individual traits (e.g. Langan, Higgins, & Scheiter, 2017; Sakschewski multi-dimensional trait space of functionally related traits (Chave et et al., 2015) which explore and map the mechanisms between FD al., 2009; Díaz et al., 2016; Mason, Mouillot, Lee, & Wilson, 2005; and ecosystem functions from local to regional scales. The interplay THONICKE et al. | 1071 between flexible morphological and physiological traits within plant 2 | M ATE R I A L S A N D M E TH O DS functional types (PFT, Lavorel et al., 2007; Prentice et al., 2007) in combination with the physiology and biogeochemistry of a DGVM We connected the leaf and stem economics approach as imple- allows analysing the effect of community assembly on ecosystem mented in LPJmL-FIT (Sakschewski et al., 2015) with a phenology functions, for example, productivity and carbon storage, in forest model (Forkel et al. (2014) to simulate potential natural vegeta- ecosystems. Because these flexible-trait DGVMs vary plant traits for tion in Europe under current climatic conditions. Herbaceous PFTs individual trees that belong to a specific PFT, inter-PFT as well as (C 3 and C4 grasses) were simulated as in LPJmL (Schaphoff, von the intra-PFT trait diversity are captured which allows investigating Bloh, et al., 2018). To allow for environmental and competitive fil- effects of niche complementarity along climatic gradients. tering to take full effect within and across PFTs, we removed the The overall aim of this study is to investigate the interaction bioclimatic limits that are used in most DGVMs to emulate bio- among climate, ecosystem functions and pattern of FD of forest geographic limitations of PFT occurrence (Schaphoff, von Bloh, et ecosystems. Therefore, we adapted the DGVM Lund-Potsdam- al., 2018; Sitch et al., 2008). Phenology, leaf-economics (LES) and Jena managed Land of flexible individual traits (LPJmL-FIT) to PFTs stem-economics (SES) traits (cf. Chave et al., 2009; Wright et al., growing in strongly seasonal European climatic conditions, while the 2004) were assigned to each individual tree sapling at establish- model previously has been successfully applied to tropical rainfor- ment allowing any trait combination (everything-is-everywhere ests (Sakschewski et al., 2015). We break down the overall aim of the approach, see Figure 1 for trait-selection algorithm) whose com- study into the following research questions: petitiveness in a given climate then determines its survival and growth. The implemented LES traits include Specific Leaf Area 1. What is the role of environmental and competitive filtering on (SLA, leaf area per unit leaf mass, mm2 mg−1), Leaf Longevity (LL, trait distributions and productivity of Mediterranean, temperate average life span of leaves, in months), leaf nitrogen content (Nleaf, and boreal natural forests? leaf nitrogen content per leaf area, mg/g) and Vcmax (maximum 2. How does FD emerge from climate and plant competition at the carboxylation rate of the RUBISCO enzyme per leaf area at 25°C, local and pan-European scale? μmol CO2 m−2 s−1). While the LES traits influence vegetation pro- 3. How does FD vary between and across European natural forests? ductivity, the SES trait Wood density (WD, wood dry mass per unit of green volume, g/cm3) is related to biomass and tree mortality We focus on European natural forests ranging from broad- (for more details, see Sakschewski et al., 2015). Physiology, growth leaved evergreen vegetation in the Mediterranean basin, to tem- and mortality of trees within the forest patch were as described in perate forests and to boreal forests in northern Europe. Here, Sakschewski et al. (2015). Applying this simulation framework to natural forests are defined as potential natural forests whose current European climate, environmental and competitive filtering compositions and ecosystem functions (see Geller et al., 2017; result in site-specific trait distributions, productivity, biomass and Hooper et al., 2005 for definition) results from climate and soil in tree height. conditions. Forests and other wooded land are usually defined fol- LPJmL-FIT combines flexible individual traits with gap dynam- lowing the percentage of woody cover (FAO, 2018). However, we ics and plant physiology, hydrology and biogeochemistry. Being denote forests hereafter as vegetation with a significant amount structured into vertical leaf layers every two meters, trees com- of biomass (>50 gC/m2), at least 5% coverage of woody PFTs and a pete for light and water as they grow in size. The trait combination minimum mean tree height of 5 m. of each tree determines its competitive strength under given cli- To address these research questions, we check the validity of the mate conditions at a given site, where several plant strategies can adapted LPJmL-FIT model by assessing to what degree the model re- co-exist and form diverse communities in forest ecosystems. The constructs observed (a) productivity and biomass, (b) trait distribution suitability of the trait combinations to the local climate conditions and (c) distribution of PFTs as a result of environmental and compet- determines trees’ competitiveness, and thus which tree strategy itive filtering. We quantify for selected sites and on a Pan-European is dominant from wet to dry tropical conditions (Sakschewski et scale FR, FDv and FE of simulated physiological and morphological al., 2015). In general, the model approach of LPJmL-FIT allows all traits (cf. Schneider et al., 2017; Villeger et al., 2008). We expect FD tree strategies to establish everywhere at any time (“everything to be relatively high in natural forests as was found in, for example, is everywhere”). Trees that underperform due to trait combina- Schneider et al. (2017) and climatic stressors considered in LPJmL-FIT tions less suitable for local climate and/or competitive conditions to be a strong environmental filter (Bernard-Verdier et al., 2012). have a relatively high mortality probability (see Sakschewski et al., Our analysis of FD in natural forests can provide a reference 2015 and 2.1 below). Under stable climate conditions, establish- state for restoring highly managed or degraded ecosystems and in- ment and mortality lead to an equilibrium of tree abundance and crease their diversity and stability in face of climate change (Mori, community composition in which only those trees occur whose Lertzman, & Gustafsson, 2017). We focus here on the interaction functional traits are suited to the local environmental conditions of physiological and morphological plant traits with community as- (climate and soil). sembly processes at the local and Pan-European scale to understand In the following, we describe the adjustments required to apply how FD emerges from those interactions. the LPJmL-FIT model to European natural forests. 1072 | THONICKE et al. F I G U R E 1 Conceptual scheme illustrating how flexible individual parameter are assigned to tree individuals belonging to a specific woody PFT (Steps 2 and 3), including LES (leaf-) and SES (stem economics) plant traits in LPJmL-FIT (SLA – Specific leaf area; WD – Wood density). Parameter sets consist of Phenology (temperature (T), light and water), LES and SES traits (WD) and are assigned to individual tree saplings at model initialization and gap opening (step 4). Climate, soil conditions (environmental filtering) and competition between individual trees within and across PFTs (competitive filtering) result in aggregated grid-cell trait distributions within patches (top-right figure in output panel). They are visualized as trait maps and are used to quantify diversity indices. PFT: BL-S – Broad-leaved summer green tree, BL-E – Broad-leaved evergreen tree, B-NL – Boreal needle-leaved evergreen tree, T-NL – Temperate needle-leaved evergreen tree. LPJmL-FIT, Lund-Potsdam-Jena managed Land of flexible individual traits; PFT, plant functional type 2.1 | Adjusted LPJmL-FIT model boreal forests under the influence of light, water and tempera- ture stress, we implemented the phenology model of Forkel et Lund-Potsdam-Jena managed Land of flexible individual traits al. (2014) into LPJmL-FIT. This model couples the phenological was originally developed for Amazonian rainforests (Sakschewski status of a tree, which ranges between 0 (complete senescence) et al., 2015). Adapting it to climatic conditions which are strongly and 1 (fully leaved), to the local climate. The actual value of the seasonal and with a high intra-annual variability requires a num- phenology status is determined by the product of four phenology ber of adaptations and implementations for multiple, co-oc- functions, which depend on a set of PFT-specific parameters and curring PFTs to describe Mediterranean, temperate and boreal the daily temperature, water stress and radiation. We calibrated natural forests in Europe. The LES and SES approach, as imple- the phenology parameter (Table S1) to yield a best possible PFT mented in LPJmL-FIT, has therefore been adapted to accommo- distribution that matches the spatial distribution of European date four tree PFTs, namely, “broad-leaved summer green” (BL-S), natural vegetation from Bohn et al. (2007). In addition, leaf se- “broad-leaved evergreen” (BL-E), “Temperate needle-leaved ev- nescence now occurs immediately if the phenological status of ergreen tree” (T-NL) and “Boreal needle-leaved evergreen tree” a tree drops below 0.2, forcing a tree to rebuild the complete (B-NL). Each PFT is based on an earlier implementation of these canopy leaf area in the next simulation year. Given the climate PFTs in LPJmL-4 (Schaphoff, von Bloh, et al., 2018) which has influence on leaf phenology changes along continental climate been extensively evaluated (Schaphoff, Forkel, et al., 2018). To gradients, a continuous spectrum from summer to winter decidu- account for the phenology of the Mediterranean, temperate and ousness for BL-S trees is captured in the model. THONICKE et al. | 1073 The algorithm to combine the PFT parameter set of the LES- and at European sites still containing natural forests (Hainich National Park SES-related traits (WD) with the phenology parameter is implemented (NP) and Bialowieza NP) because the model simulates natural vegetation as follows: When a new sapling establishes, the selection of the trait only. Following this calibration, we set the parameter 𝛼1 to −4.5 and 𝛼2 combination occurs in three steps. First, the PFT type is randomly cho- to −2.66 for the broad-leaved trees, and 𝛼1 = −2.66 and 𝛼2 = 0.255 for sen out of the four possible PFTs independent of its climate suitability needle-leaved trees. The term mortWD is used as the maximum of the or the present PFT composition (Figure 1). Second, the selected PFT growth-efficiency mortality in LPJmL-FIT, meaning that trees with a low type defines the set of phenology parameters that are assigned to the growth efficiency resulting from low productivity under unfavourable tree sapling (see Table S1), and it defines the SLA range from which the climate conditions have a higher mortality risk (Sakschewski et al., 2015). SLA value is drawn from a uniform distribution. The SLA range differs In addition, we reduced the tree allometry parameter kep to 1.5 for BL-E, T-NL, B-NL and BL-S (see Table S2). The functionally related (1.6 in Schaphoff, von Bloh, et al., 2018) for needle-leaved trees to plant traits leaf longevity (LL), leaf N and Vcmax (maximum carboxyl- simulate realistic tree shapes and growth pattern of needle- in com- ation rate of Rubisco per leaf area) are assigned in a third step to the parison to broad-leaved saplings: sapling following the LES approach as in Sakschewski et al. (2015), see Figure 1. Both, SLA and WD are drawn from a continuous uniform dis- CA ∼ Dkep (3) tribution for each individual tree for which the PFT-specific range is derived from the TRY v4.0 database (see https​://www.try-db.org/ and where kep mediates between crown area (CA) and stem diameter (D). Kattge et al., 2011) considering only sites located inside our study area. The lower kep, the CA for needle-leaved trees is reduced which then Hence, the parameter set of each newly established sapling in the for- affects LAI: est patch contains trait values drawn from the LES and SES as well as the phenology parameters. This approach is applied to the 20 new sap- Cleaf ⋅ SLA LAI = , (4) lings that are established per patch each year. Grasses are integrated CA in the model as homogeneous layers using the remaining radiation at the bottom of every forest patch for photosynthesis. C3 and C4 grasses where Cleaf is whole plant carbon investment in leaves (kg per tree) (cf. compete with trees for water only. The establishment rate of grass is Sitch et al., 2003). Because of the lower SLA, needle-leaved trees then anti-proportional to the actual tree cover and the mortality rate de- have to invest more leaf carbon in their first years to reach the same pends on carbon balance by the end of every year (see Schaphoff, von LAI compared with broad-leaved trees. Bloh, et al., 2018 for details). Grasses are assigned phenology parame- Fire is an important natural disturbance in European forest ecosys- ters (see Table S1), but are simulated without individual trait flexibility tems (Naveh, 1990; Tinner et al., 1999). We applied the simple Glob- following the LPJmL4 modelling approach (Schaphoff, von Bloh, et al., FIRM model (Thonicke, Venevsky, Sitch, & Cramer, 2001) as embedded 2018). Therefore, in this study model simulations will be used to evalu- in LPJmL4. Here, fire probability depends on soil moisture in the top soil ate trait distributions and FD regarding trees only. layer and a fuel load threshold modelled by an exponential probability Different from the LPJmL-FIT version of Sakschewski et al. function calculated by the end of every year. We apply this probability (2015), the trade-off between SLA and LL has been adopted from to each patch separately. If a patch is burnt, every tree is ignited and LPJmL4 (Schaphoff, von Bloh, et al., 2018): survives with a PFT-specific fire resistance probability (cf. Thonicke et al., 2001) which is the same for all individuals that belong to the same PFT. In the patch burnt, all litter carbon is combusted completely. ( ) SLA ⋅DM (1) 𝛽0 −log10 ∝ c LL = 10 𝛽1 In boreal forests, evergreen trees exhibit water stress in spring, when relatively high air temperature increases evaporative demand, where DMc denotes the dry matter carbon content of leaves while the soil is still frozen and thus limits the water availability in the DMc = 0.4763 and the parameter ∝= 2 ⋅ 10−4. The parameter 𝛽1 is set soil. This water stress forces evergreen trees to shed their needles and ( ) to 0.4 and 𝛽0 to 2.2 for broad-leaved PFTs, while 𝛽0 = 2.08 for both, making them less competitive against BL-S trees which might have T-NL and B-NL. All parameters in Equation (1) were obtained from their bud burst later in the year. To correctly balance competition be- Kattge et al. (2011). The parameters 𝛽1 and ∝ influence the steepness tween B-NL and BL-S in boreal forests, we increased the root-distribu- of the SLA-LL relation, whereas 𝛽1 alters the offset. tion factor 𝛽root for B-NL trees from 0.943 (Schaphoff, von Bloh, et al., The mortality mortWD is coupled to its wood density WD by 2018) to 0.965, which allows B-NL to reach deeper, non-frozen layers using the equation from King, Davies, Tan, and Noor (2006): during spring thereby preventing them from leaf senescence. mortWD = 10𝛼1 +𝛼2 ∕WD (2) 2.2 | Model input data, simulation protocol and and assigned to each tree individual at establishment (Sakschewski et al., validation sites 2015). Because no general mortality–WD relationship for tree species of temperate forests is currently available in the literature, we calibrated To post-process and analyse our simulation data, we used R3.6.0 𝛼1 and 𝛼2 to the locally observed biomass. Calibration was carried out (R-Core-Team, 2019) and MATLAB 2019 (MATLAB, 2019). All 1074 | THONICKE et al. R-packages used in this study including their citation are listed in remote-sensing data (MOD17A2H) for the years 2004–2013 Table S6. (Running, 2015) at six sites (see Table S1) and respective flux tower measurements from the Euroflux network for the Laegeren (CH- Lae, D'Odorico, 2014; Paul-Limoges, 2018) and Hainich NP (DE- 2.2.1 | Climate and soil data input Hai) (for general information, see Papale et al., 2006; Reichstein et al., 2005). Simulated maps of vegetation height and biomass were Lund-Potsdam-Jena managed Land of flexible individual traits uses evaluated against remotely sensed products (Lefsky, 2007; Thurner air temperature [°C], precipitation [mm/d] and radiation (short-wave et al., 2014). Details on the validation of Gross Primary Productivity down and long-wave net radiation [W/m2]) of the combined data- (GPP), vegetation height and biomass are described in the Data S1. set of the WATCH (Weedon et al., 2011) and WFDEI (Weedon et Simulated plant trait distributions were compared against observed al., 2014) datasets at daily resolution on a 0.5° x 0.5° longitude– plant trait data from the TRY database (Kattge et al., 2011), see Data latitudinal grid. This climate data set is based on the reanalysis of S1 for methods and data origin. ERA-Interim, where precipitation was bias corrected using the Global Precipitation Climatology Centre data set (GPCC, Schneider et al., 2011). The climate data range from 1901 to 2013 with WFD 2.3 | Computation of FD indices (WATCH Forcing Data) covering 1901 to 1978 and WFDEI-GPCC (WATCH-Forcing-Data-ERA-Interim) is used from 1979 onwards. We quantified three complementary indices on multidimensional The atmospheric CO2 concentration is held constant at 296 ppm traits to describe FD, namely, FR, FDv and FE, following Villeger et over the whole time period. al. (2008) and Schneider et al. (2017), where each point represents Soil texture is needed as model input and was taken from the one tree individual (higher than 2 m) with its unique trait com- Harmonized World Soil Database version 1.2 (Nachtergaele et al., bination. In this study, this multidimensional trait space is based 2009). The soil depth was kept constant at 2m for all grid cells. Our on SLA, LL, WD and tree height. While SLA, LL and WD influ- simulation domain covers Europe from 11°E to 36°W and 29.5°N to ence productivity and biomass (Reich, 2014) and therefore point 62°N. to competitive exclusion, tree height is regarded to describe niche differentiation (Garnier et al., 2016). We calculated the FD indi- ces across and within PFTs to capture assembly processes across 2.2.2 | Simulation protocol meta-communities. Functional richness describes the extent of the occupied trait Model simulation starts from bare-ground and simulates a spin-up space and is calculated by the convex hull volume including all period of 500 years by recycling the first 30 years of the climate points in that trait space, which is normalized by the maximum data set (1901–1930) to bring natural vegetation composition (here possible trait volume. However, it implies that FR reacts strongly individual trees with their individual trait combinations) and all living to outliers. FDv describes how far environmental niches are sepa- and dead carbon fluxes into equilibrium with the spin-up climate. We rated and indicate the intensity of competitive interactions, where then performed a transient run simulating potential natural vegeta- FDv = 0 indicates convergent trait distribution due to strong en- tion until the end of 2013, that is, without land use. For the European vironmental and competitive filtering (Mason et al., 2005; Villeger simulation domain, 1,000 forest patches being equivalent to 10 ha et al., 2008). To measure FDv in a multidimensional trait space, a of forest area are simulated in each grid cell where all patches re- sphere with radius dG centred in the trait cloud is calculated. FDv ceive the same climate data and the same soil data as model input. then quantifies how points (trait combinations of trees) scatter rel- Respective model output is then aggregated over all simulated ative to the surface of the sphere (see eqs. 5–7 in Data S1 and patches within a grid cell. Figure S1). If all points are located on the sphere, FDv becomes To allow for a detailed analysis of plant-trait distribution, FD unity independent on the respective distribution on the sphere. and productivity, we chose six different sites across Europe with The more the index decreases, the wider the points are spread near-natural forest stands and which cover a broad range of climates around (inside and outside) the sphere. FDv therefore quantifies (Table S4). Site-specific simulations follow the same protocol as de- how the occupied niches are separated. FE describes how regu- scribed above and were performed with 2,500 patches at each site larly points are distributed in the trait space, that is, how efficient to ensure a higher spatial coverage. available resources are used through niche occupation. FE is based on the minimum spanning tree (MST) linking all points in the trait space in such a way that the sum of all branches becomes minimal 2.2.3 | Model validation (see eqs. 11–13 in Data S1 and Figure S1). Therefore, FE increases when a) the points are evenly distributed, that is, having equal Simulated seasonal and intra-annual GPP is validated against ob- branch length; or b) the trait combinations of the trees are equi- served and remotely sensed GPP data at six sites (Table S3) cov- distant in the trait space (Villeger et al., 2008). Further details on ering a climatic gradient (Table S4). We used monthly MODIS FE, FDv and FR quantification are provided in Data S1. THONICKE et al. | 1075 Each index was computed by considering all trees as evenly distributions therefore emerge from the LPJmL-FIT modelling frame- weighted. Before calculation, traits were normalized to their min- work (Figure 1). We compare simulated SLA and WD against TRY imum and maximum observable values in the TRY sites in Europe, observations for BL-S, BL-E trees and needle-leaved evergreen, that thus ranging between 0 and 1. For each grid cell we checked that no is, B-NL and T-NL trees (Figure 2). Simulated mean trait distribu- dimensional reduction was required by using the function “dbFD” tions match the TRY observation reasonably well for both, SLA and of the R-package “FD” performing a principal coordinates analysis WD, for BL-S, BL-E and the two needle-leaved PFTs (dashed lines in (Laliberté et al., 2014). Due to constrained computation capacity, Figure 2). Simulated ranges of SLA, however, are smaller than the we calculated all FD indices separately in groups of 50 patches in original trait range (see Table S2) and smaller than observed SLA each grid cell (for which 1,000 patches were simulated in total) and (Figure 2). LPJmL-FIT simulates a mean SLA of 16.11 with a stand- aggregated them to the grid level by using the arithmetic mean. To ard deviation of ±1.25 mm2/mg for BL-S compared to 14.95 ± 6.09 visualize the stochastic uncertainty of the model, we calculated the documented in TRY. The simulated range for SLA is also smaller for coefficient of variation (COV) of each index in a grid cell out of the BL-E (LPJmL-FIT: 9.31 ± 1.13 mm2/mg; TRY: 6.95 ± 2.72 mm2/mg) groups of 50 patches (n = 20). As all of these groups in a cell received and the needle-leaved evergreen (LPJmL-FIT: 6.11 ± 0.95 mm2/mg; the same climate data, the COV can be seen as a measure for the TRY: 6.95 ± 2.72 mm2/mg). Simulated ranges for WD are quite close stochastic uncertainty of the model. to observed ranges in TRY for BL-S, slightly smaller for BL-E, but broader for the needle-leaved PFTs (Figure 2, bottom row). The combined effects of environmental filtering and plant com- 3 | R E S U LT S petition for light and water in LPJmL-FIT also result in reasonable seasonal and interannual productivity (GPP) as observed on six se- 3.1 | Climate influence on trait distribution, lected sites. The model quality is shown by a high Pearson's R (≥0.89) productivity and tree height and low NMSE (0.03–0.29; see Table S3, Figure S2 and Data S1 for details on the evaluation methods and results, sites are described Environmental and competitive filtering allows those trees to es- in Table S4). The simulated competition between tree individuals tablish and survive whose trait combinations are suitable for local results in closed forest cover and corresponding high biomass stor- climate conditions in the forest patches simulated by LPJmL-FIT. age in temperate and boreal forests (Figure S3). We found simulated Unsuitable or less suitable trait combinations lead to a low growth vegetation height and biomass to compare well against remotely efficiency and are therefore outcompeted. The combination of cli- sensed observations and local in situ data, although the comparison mate suitability and competiveness has the effect that the continu- of natural forest and actual vegetation is limited as remote-sensing ous, uniform distribution with which the model is initialized results products detect properties of actual vegetation cover which are in- in a normal trait distribution at the local scale. The resulting trait fluenced by current land use and forest management (see Data S1 F I G U R E 2 Trait distribution of Specific Leaf Area (top panel) and Wood Density (bottom panel) for broad-leaved summer green PFT (BL- S, left column), broad-leaved evergreen PFT (BL-E, central column) and both, boreal and temperate, needle-leaved evergreen PFTs (NL, right column) as simulated by LPJmL-FIT (light red) and observed by TRY database entries (light cyan). Dashed vertical lines show the simulated (light red) and measured (light cyan) average trait values for each PFT. LPJmL-FIT, Lund-Potsdam-Jena managed Land of flexible individual traits. PFT, plant functional type 1076 | THONICKE et al. for details on the evaluation method and results). Centuries of land Abiotic conditions, here aggregated to MAT [°C] and MAP clearing, agricultural and forestry have greatly changed land cover [mm], are strongly linked to tree establishment (Figure 3). The cli- and reduced natural forests to few remaining small areas (Ellis et al., mate space occupied by trees across all PFTs converges towards 2013), which further complicates the evaluation of simulated poten- higher MAT and lower MAP. Intra-PFT variation in SLA values tial natural vegetation. decreases with warmer and drier conditions. Most SLA variation, however, happens between PFTs that occupy different parts of the climate space (Figure 3a). The SLA of BL-S varies from values 3.2 | Climate influence on trait distributions around 13 mm2/mg under warmer and drier climate conditions to within and across PFTs >20 mm2/mg in colder climate conditions (<7°C) with a wide range in precipitation (500 to >2,000 mm MAP). BL-S co-exists with B-NL Spatial distribution of simulated fractional cover of each PFT result in cold/wet climate conditions, with BL-E in warmer and increas- from the PFT-specific phenology and functional trait combinations, ingly drier climate conditions (>10°C MAT and 300–2,000 mm which determine the suitability of each parameter set to the climate MAP). Although many trait combinations are possible under the and soil conditions in a given grid cell (cf. Figure 1). Tree individu- “everything-is-everywhere” approach of LPJmL-FIT, BL-S with SLA als with trait combinations adapted to current climate are the most values lower than 13 mm2/mg do not occur despite a possible min- competitive and most productive in sites with limited environmental imum of 7 mm2/mg (see Table S2). BL-E cover a similar tempera- stress severity (while their productivity may decrease at stressed ture range as BL-S, but occur across a wider SLA range between sites, Zhang, Niinemets, Sheffield, & Lichstein, 2018), and thus a MAT of 8–17°C, although the possible SLA range of BL-E is cover larger proportions of a given grid cell (Figure S4). The most among the smallest of all simulated PFTs. Again, this realized trait suitable combination, to shed leaves under cold and/or dry climate space separates these PFTs clearly from the two needle-leaved conditions, results in BL-S dominating central Europe, even though PFTs, T-NL and B-NL, where T-NL increasingly dominates at MAT it also occurs – albeit at much smaller fractions – in the boreal and >10°C and MAP <1,000 mm with SLA ranging from 5 to 8 mm2/ the Mediterranean forests. B-NL dominates northern Europe, and mg, whereas B-NL shows lowest SLA values at a MAT below 10°C T-NL the Mediterranean basin, where it co-occurs with BL-E. Several (see Figure 3a; Figure S5 shows PFT-specific SLA maps, Figure S8 PFTs co-occur in the Mediterranean forests (three tree PFTs and C3 maps the COV of SLA). grasses), whereas temperate forests in lowland Europe are domi- Simulated mean WD clearly separates along climate gradients nated by just one tree PFT (Figure S4). Note that in all cases, the tree across all PFTs (Figure 3b). Low WDs are simulated in cold climate individuals still vary in their trait combination within each PFT. conditions, and WD increases with increasing MAT and decreasing F I G U R E 3 Mean Specific Leaf Area (SLA) and Wood Density (WD) distribution for each PFT plotted against mean annual temperature (MAT) and mean annual precipitation (MAP). The colour of each dot represents a mean SLA (a) and WD (b) value and each panel belongs to a woody PFT. (Broad-leaved summer green: BL-S; Broad-leaved evergreen: BL-E; Boreal needle-leaved evergreen tree: B-NL; and Temperate needle-leaved evergreen tree: T-NL) Mean trait values are averaged over the 2004–2013 time period for 1,000 patches per grid cell. PFT, plant functional type THONICKE et al. | 1077 MAP. B-NL and BL-S cover the WD space between 0.5 and <0.7 g/ and Dundo (Figure 4, left column). The occupied trait space forms 3 3 cm below 10°C MAT. WD around 0.7 g/m are found in cold (<10°C the hypervolume, that is, the FR (shown in grey-blue in Figure 4, 2nd MAT) and drier climate conditions (<1,000 mm MAP), with similar WD column) for each site. The sphere around the centre of gravity (grey values found in warmer (10–12°C MAT) and wetter (>1,000 mm MAP) surface and green cross, respectively, shown in Figure 4, 3rd column) climate conditions. Highest WD values (>1.1 g/cm3) are simulated for illustrates site-specific FDv, whereas FE is quantified from the MST T-NL, BL-S and BL-E with >15°C MAT and <1,000 mm MAP (Figure 3b). (Figure 4, last column). Table S4 shows the site-specific FD indices for Higher WD allows slow plant growth and lowers tree mortality risk the six sites. (Equation (2) in 2.1) which explains why high WD values are simulated The wide bi-modal distribution of SLA between BL-S and B-NL across all PFTs under dry climate conditions (see also Figure S6 for trees in Seitseminen increases the trait space, that is, FR, whereas PFT-specific WD maps and Figure S9 for maps displaying the COV of in Laegern and Dundo simulated SLA distributions show narrower WD). bi-modal distributions or even converge (see density distribution in Figure 4, left column). Niche separation (FDv) and regularity of niche occupation (FE) are more comparable across the three sites (Table S4). 3.3 | FD emerging from climate and plant FDv is highest in Dundo because points in trait space lay closer to the competition surface of the sphere compared with Seitseminen and Laegern (notice points outside the sphere in Seitseminen and Laegern). Compared with The calculation of all diversity indices is based on the initial four-di- Laegern, we find slightly higher FDv in Seitseminen because of the di- mensional trait space out of SLA, LL, WD and tree height. However, vergent SLA distribution. Niche occupation is less regularly distributed for visualization, we remapped the trait space from four dimen- (FE) in Dundo compared with Seitseminen and Laegern because the sions to a three-dimensional trait space composed of SLA, WD and trait space of B-NL trees is less occupied in Dundo, leading to larger tree height. We plot the position of each tree individual in the trait path length in between points of this PFT (Data S1, Figure S1). This space for the climatologically different sites Seitseminen, Laegern leads to more irregular distances in between points, which lowers FE. F I G U R E 4 Distribution of simulated trees in a three-dimensional trait space (Height, WD and SLA). We plot density (far left) and functional diversity indices for 2,500 simulated patches (FR centre left; FDv centre right and FE (far right) for three sites (Seitseminen (top row), Lägern (middle row) and Dundo (bottom row)). Ranges of the SLA and WD axes correspond to the maximum trait range across all woody PFTs used in the simulation (see Table S2). Seitseminen, Laegern and Dundo represent boreal, mixed temperate and Mediterranean- type forests respectively. Each dot represents a trait combination of one tree larger than 5 m while the colour indicates its PFT type. Plant Functional Types are as followed: Broad-leaved summer green (BL-S), Broad-leaved Evergreen (BL-E), Boreal needle-leaved evergreen tree (B-NL) and Temperate needle-leaved evergreen tree (T-NL). FE, functional evenness; FR, functional richness; PFT, plant functional type; SLA, specific leaf area; WD, wood density 1078 | THONICKE et al. When calculating the FD indices for each PFT trait space separately on the British Isles. Where B-NL and BL-S occupy distant parts of the and for each site, FDv and FE are similar to the corresponding across- trait space, FR reaches its maximal values of 0.03–0.04, very much PFT values (Table S5 and Data S1 for computation of within-PFT FD alike in the Seitseminen site (compare Figure 4, top row). In contrast, indices). Within each PFT, diverse functional strategies co-exist, that lowest FR values are computed for areas where one PFT is dominant, is, the intensity of plant competition and regularity of niche occupa- especially in the lowland areas of temperate forests which are domi- tion are comparable to the one across PFTs. Specifically, FDv is high nated by BL-S trees (Figure 5a). The COV for FR is high in temperate at Seitseminen for B-NL trees compared with the overall FDv because and alpine forests, where mean FR is low because variability increases the point cloud is clearly separated along the tree-height niche axis where the mean of a variable is close to 0. On the contrary, low values and shows a wider WD distribution compared with the other sites are found in Mediterranean and boreal forests (Figure S7). (Figure 4, left column). However, intra-PFT FR is three orders of magni- Functional divergence is high in natural forests reaching val- tude lower than FR across PFTs. This is mainly caused by a much lower ues between 0.68 and 0.82 (Figure 5b). Where needle-leaved and realized trait space for LL (e.g. BL-S) and SLA (e.g. B-NL, see Figure 4), broadleaved trees co-exist (cf. Figure S4), FDv is higher, that is, in bo- whereas the intra-PFT range of tree heights and WDs is similar to that real and mountain forests, and in southern Mediterranean forests. between PFTs. In summary, environmental and competitive filtering in- Where only one PFT dominates, FDv is lower (0.7–0.73 compared fluence niche occupation in a similar way within as well as across PFTs. to >0.75), for example, in lowland temperate and Mediterranean forests. In the transition zone to boreal forests or mountain forests, trait distributions of BL-S-dominated forests further converges 3.4 | FD at the European scale (FDv ~0.68). Here, intensity of filtering increases (FDv converges) when climate conditions reduce growth efficiency of BL-S trees. At the European scale, spatial gradients in FD indices are relatively Further north, growing conditions for B-NL are more suitable, al- small and a few spatially distinct patterns stand out (Figure 5). FR in- lowing establishment of another plant strategy causing an increase creases with the number of PFTs present through which the size of the in FDv (Figure 5b). In contrast, in the southern Mediterranean for- trait space increases (Figure 5a, compare Figure S4). Higher FR is found ests, needle-leaved and broad-leaved trees are smaller (Figure S3) in mountain areas throughout the continent, in boreal forests but also and in competition with grasses for water, thus, FDv is higher again F I G U R E 5 Maps of Functional Richness (FR, panel a)), Functional Diversity (FDv, panel b)) and Functional Evenness (FE, c)). All three functional diversity indices are computed from specific leaf area, leaf longevity, wood density and tree height averaged over all individual trees, which are at least 5 m tall. Effects on FR, FDv and FE result from trait variability within and across woody PFTs present in 1,000 patches within a grid cell. PFT, plant functional type THONICKE et al. | 1079 (~0.76, see Figure 5b). COV of FDv is generally very low with max- of trees surviving and growing in a forest patch. Additional pro- imum values found in some mountainous and boreal transitional cesses not yet explicitly captured by LPJmL-FIT might explain the areas (Figure S7). remaining trait variability and include: (a) dispersal and adaptive Functional Evenness follows a different spatial pattern across responses (incl. phenotypic plasticity of traits), (b) nutrient avail- Europe than the other two indices (Figure 5c). Lower FE values ability (e.g. nitrogen limitation), (c) variable rooting strategies and (~0.72) are found in lowland temperate forest in central and east- disturbances other than fire (Van Bodegom et al., 2012; Douma ern Europe as well as in the transition zone to the boreal forest in et al., 2012) and (d) trade-offs between different trait combina- north-eastern Europe, where BL-S dominate but still grow in com- tions and species capacities to tolerate multiple environmental petition with B-NL trees (Figure S3). FE increases where only one stresses simultaneously (Laanisto & Niinemets, 2015; Niinemets NL-PFT dominates (e.g. British Isles, Norway) and is high in wet & Valladares, 2006). While the influence of adaptive responses on Atlantic and in southern Europe, where high environmental stress trait distribution is perhaps difficult to measure, including nutri- results in efficient resource use. Here, COV of FE corresponds to the ent availability would introduce another niche axis and perhaps spatial gradient of mean FE, maximum COV is found where mean FE differentiate the trade-off with SLA in the LES better. Variable is lower (Figure S7). rooting strategies could alleviate water stress in seasonally dry environments and increase growth efficiency, hence tree height of Mediterranean forests (cf. Figure S3). We have included fire 4 | D I S CU S S I O N disturbance in LPJmL-FIT, which simulates high fire activity in the Mediterranean forests, less in boreal and low fire activity in tem- 4.1 | Effect of environmental and competitive perate forests, but other disturbance agents such as windthrow filtering on trait distribution and productivity or frost could also further diversify plant strategies, and thus in- crease the width of the simulated trait distribution. Implementing The adapted LPJmL-FIT model is capable of reproducing observed these disturbances into LPJmL-FIT would be expected to decrease GPP with a small modelling error and high correlation with ob- the competitiveness of the BL-S trees in boreal and mountainous served data (see Data S1). Biomass and plant height follow the forests, and thus further restrict their spatial extent. spatial distribution of previous publications (Thurner et al., 2014; Specific model functions determine the model outcome. The Healey et al., 2015), although the comparison is limited by the long- phenology parameter (Table S1) influences GPP, growth efficiency term land use history in Europe which restricts the comparability and thus the performance of a tree in the forest patch. In exchange between simulated and observed data and may explain the discrep- with climate conditions, the phenology parameters determine the ancy found. Simulated SLA and WD are in close agreement with spatial distribution of the PFTs, including grasses. The trade-off observed TRY data. These results show that the new version of between SLA and LL defines the carbon cost for leaves (or nee- LPJmL-FIT reproduces the spatial PFT and local trait distributions dles) under specific climate conditions which then determine the as well as the productivity and biomass of European natural for- length of the growing season, and thus vegetation productivity. ests. Even though a preference for measuring short-lived, broad- With more SLA and LL measurements available for temperate for- leaved species can influence measured trait data which could ests, this relationship could be adjusted and contribute to further overestimate, for example, SLA (Sandel et al., 2015), such bias is reduce modelling error of the simulated trait distribution (Figure 2). found to be small for European trees and should have little influ- The WD-dependent mortality function (see Equation (2)) influences ence on our results because LPJmL-FIT clearly separates between simulated biomass, therefore more data on WD and tree mortality broad-leaved and needle-leaved trees and focuses on trees only. would help to improve parameterization of this function. Further The model is capable of simulating potential natural vegetation reductions in the model error of simulated WD distribution, tree without prescribing the spatial extent of PFTs via bioclimatic limits height and biomass could be expected. and without prescribing the functional traits in question (e.g. SLA, Plant-trait validation would profit from more and better resolved LL, WD). Therefore, it enables to investigate the interaction be- plant trait measurements in natural forests in Europe. In this study, tween environmental and competitive filtering in European natural we had to aggregate several TRY sites and also to merge T-NL and forests in an unprecedented manner. B-NL to validate simulated trait distributions. Ideally, trait measure- Simulated ranges for SLA (BL-S and BL-E) and WD (BL-E and ments at FLUXNET sites, where we can evaluate both productiv- needle-leaved trees) are smaller than observed (Figure 2 and Table ity and plant traits, would allow high resolution and in-depth model S2). Heat and drought stress as well as light availability (seasonal, evaluation. With these consistently measured data available, model vertical canopy structure) influence growing conditions of trees. evaluation could be extended to water fluxes and biomass as well These climate factors act as environmental filter favouring plant as stand structure and their related plant traits. Alternatively, re- strategies that are adapted to the local climate. Tree growth effi- motely sensed traits which are already available at the Laegern site ciency influences tree competitiveness and its capacity to reach (cf. Schneider et al., 2017) or emerging from global remote-sens- the top layer in the canopy and gain most light in temperate and ing missions (cf. Jetz et al., 2016; Ma et al., 2019) could be used in boreal forests. Therefore, competitive filtering is the second driver the future. However, cross-validation experiments are required to 1080 | THONICKE et al. compare spatially continuous traits observed from remote sensing PFTs are distant and thus cover a large volume of the trait space. In to traits simulated by DGVMs using flexible-individual traits, but ap- high-elevation areas and the boreal zone, BL-S adapt to high SLA proaches are underway (cf. Garonna et al., 2018; Lausch et al., 2018). values and lower LL because of shorter vegetation periods lead- Simulated trait maps for SLA and WD can be validated against inter- ing to a wider functional gap in these two traits between BL-S and polated observed trait maps (Butler et al., 2017; Šímová et al., 2018), B-NL trees. FR is small due to small trait variation, where only BL-S even though observed trait distributions include those from highly trees dominate the patches of temperate forests (Figure 5). In the managed forests. Interpolated trait maps, which combine plant trait Mediterranean region, FR is comparable to temperate forests be- data with actual species’ presence data, can lead to a potential bias cause lower tree heights limit trait space and counteract the effect in data-model comparison because species are planted in managed of diverse leaf strategies. In the southern Mediterranean, where tree forests. Generating trait maps which use potential natural species growth reaches its limits, FR decreases because few niches can be distribution would improve such data-model comparison because occupied under the dry climate conditions. they would be closer to the spatial trait distribution that would result Functional divergence describes the diversity of co-existing from environmental and competitive filtering, for example, at local functional strategies. A high FDv illustrates a high degree of niche scales (Benavides, Scherer-Lorenzen, & Valladares, 2019). Model- differentiation, thus high competitive exclusion or intense competi- data comparison could be improved with a) more trait data being tive interaction (Garnier et al., 2016; Mason et al., 2005; Villeger et extracted in natural vegetation and extrapolated using species dis- al., 2008). The larger the influence of environmental filtering through tribution data to the continental scale; and b) allowing for a better abiotic stress, the more the trait distribution converges (Bernard- overlap of observed versus simulated plant traits such as leaf N (ob- Verdier et al., 2012). Maps of FDv (Figure 5) therefore illustrate servation available but simulated map missing) or WD (observation the relative strength of environmental versus competitive filtering missing, but simulated map available). on trait composition. Where climate allows needle- and broad- leaved trees to co-exist (in mountains, boreal forests, wet temper- ate and Mediterranean forests), FDv is high because of the distant 4.2 | Community-assembly effects on traits and FD SLA distribution of BL-S and B-NL. In Mediterranean and temper- across Europe ate climates, SLA distributions converge more, and FDv declines (Figure 5b). Interestingly, FDv declines further in the transition from Community assembly at a site results from dispersal or migration, temperate forests to boreal and to mountain forests respectively. environmental filtering of trait combinations adapted to local climate Being dominated by BL-S, the performance of these trees declines and finally plant competition (Bernard-Verdier et al., 2012). In the as climate gets colder, thus the abiotic stress increases which leads adapted LPJmL-FIT model, trait combinations are drawn from ob- to a less divergent distribution. Further north (or at higher elevation) served ranges for SLA and WD. It principally means that any con- the growing conditions for B-NL improve, thus their dominance in- ceivable trait combination can establish everywhere at any time creases and increases FDv again. Physiological traits (SLA, WD and for every PFT without explicitly considering trait inheritance or LL) converge at the transition from temperate to boreal forests, indi- seed dispersal. Currently, surviving trait combinations in the model cating competitive exclusion. In contrast, in southern Mediterranean represent those trees whose traits are best suited to local climate forests, increasing drought stress reduces tree height and favours a and the competitive conditions in the established tree community. wide SLA range and higher WD (Figure 3, Figure S5), and thus the Therefore, environmental filtering interacts with competitive filter- trait distributions diverge again, promoting co-occurring alternative ing and reproduces observed trait distributions and productivity of strategies (Figure 5b, Figure S3). We find that such bi-modal changes natural forests. Scaling up to the whole study region, SLA and WD in FDv reflect shifts between competitive exclusion, linked to physi- gradients emerge along climatic gradients and reflect variability in ological traits (SLA, WD, and LL), and niche differentiation, linked to trait ranges. morphological traits (tree height). Simulated SLA cover a wide temperature and precipitation range, Functional evenness quantifies the regularity of the distribution separated by SLA ranges as observed for the four PFTs in Europe. in the trait space, that is, niche occupation (cf. Mason et al., 2005; Where cold temperatures and light increasingly limit productivity, Villeger et al., 2008) and could be interpreted as a lower utilization that is, in the boreal zone, simulated SLA for BL-S increases, meaning of resources due to the more irregular occupation of the trait space, trees with extremely high SLA and short LL survive cold winters and that is, environmental niches (cf. Mason et al., 2005). However, grow in short summers. The shift towards higher WD in southern there are few studies which have investigated the changes in FE Europe indicates a local adaptation to seasonally dry Mediterranean- along climatic gradients. Recent studies focussed on changes in FE type climate. In general, the increasingly dry climate filters higher along disturbance gradients (e.g. Mouillot, Graham, Villeger, Mason, WD for BL-E and T-NL which indicates better adaptation to drought. & Bellwood, 2013) and at the local scale (Pakeman, 2011). We find The identified spatial FD patterns reflect the combined effect high FE in areas where dry and cold climate conditions limit tree of environmental and competitive filtering. The influence of climate growth and under wet (Atlantic) conditions. Increasing disturbances on surviving plant strategies and their coexistence are reflected in (fire) and climatic stress (drought) increase resource-use efficiency the occupied trait space. FR is high where trait clusters of different as suggested by Pakeman (2011). In highly stressed environments, THONICKE et al. | 1081 no strategy can dominate, thus functional traits are more evenly specialized flexible growing strategy, very growing strategy, specialist in warm, dry TA B L E 1 Overview of the woody PFT characteristics towards phenology, fire sensitivity, leaf and growth strategy, and environment. Plant Functional Types are as follows: Broad-leaved distributed. Low FE is found in mixed temperate forests and in the strategy, specialist in cold climates with Conservative leaf strategy, conservative Conservative leaf strategy, fast growing growing strategy competitive in warm Less flexible leaf strategy conservative transition from temperate to boreal forests, where BL-S and B-NL Very flexible leaf strategy, but less occupy distant parts of the trait space. At the Pan-European scale, competitive in temperate zone and fire-prone environments FE does not follow a unique climatic gradient which could be inter- preted as habitat filtering as suggested by Pakeman (2011). In our short growing season study, it is rather a combination of environmental and competitive temperate zones filtering within and across PFTs. We computed FR, FDv and FE from four traits following the ap- Environment proach of Villeger et al. (2008) as implemented for the Laegern site in Schneider et al. (2017). The FD indices change with the traits con- sidered and the dimensionality of the trait space. At the Laegern site, summer green (BL-S), Broad-leaved Evergreen (BL-E), Boreal needle-leaved evergreen tree (B-NL) and Temperate needle-leaved evergreen tree (T-NL) we compute FR values one order of magnitude lower than published Flexible WD, low WD in north, in Schneider et al. (2017) because we calculate FR for a four- instead Fast growth, due to lower WD High WD in dry environment, High WD in dry environment climate and thus lower WD fast-growing under colder of a three-dimensional trait space. In contrast, we derived compara- ble FDv and slightly lower FE values for the Laegern site (Schneider lower WD under wet et al., 2017). A direct comparison between remotely sensed and sim- high WD in south Growth strategy ulated FD based on three traits is of interest and relies on directly conditions comparable trait spaces and indices. Quantification of FD indices strongly depend on the type of traits (physiological vs. morpholog- ical) and the context of their interpretation, for example, niche dif- ferentiation or ecosystem functions to which the considered traits Mostly thin and short-lived leaves Leaf strategy (SLA range, SLA-LL Expensive and long-lived needles Expensive and long-lived needles relate. Comparing site-specific FD indices based on similar traits where also trait observations are available would be a great step intermediate leaf longevity forward to evaluate site-specific FD and its changes along climatic Moderate leaf cost and gradients at the regional as well as continental scale. Model-data Abbreviations: LL, Leaf Longevity; PFT, plant functional type; SLA, specific leaf area; WD, wood density. comparison can help to further investigate the importance of spatial scale or abiotic gradients at the landscape scale, where DGVMs with flexible individual traits such as LPJmL-FIT can be used to test and relation) explore respective FD hypotheses to advance biodiversity–ecosys- tem functions theory (Hisano, Searle, & Chen, 2018). Trees occupy trait ranges within each PFT allowing different Medium PFTs to co-exist as a result of FD, tree demography and disturbance High Low Low Fire in all biomes. Site-specific water availability and temperature have been shown to influence functional composition in temperate for- ests emphasizing the importance to investigate FD along climate Drought gradients (Zhang et al., 2018). High High Low Low We find high FD in European natural forests not only for temper- ate forests as found by Liebergesell et al. (2016) but also for boreal Warm temperature and Mediterranean forests. The simulated community assembly of European biomes can be explained by the interplay between envi- ronmental and competitive filtering: Medium High High Low 1. In the boreal zone and in high-mountain areas (e.g. Carpathian Mountains, Alps or the Pyrenees) the growing season is short, and needle-leaved trees (low SLA), which keep their leaves longer Cold temperature (high LL) but grow relatively fast due to low WD, are the dom- Resistance to inant tree growth strategy mostly as a result of environmental filtering. B-NL trees have high resistance to cold temperatures High Low Low Low and their evergreen strategy allows them to fully exploit the short growing season (see Table 1). In contrast, broad-leaved B-NL trees (BL-S), which are more prone to cold temperatures and T-NL BL-E BL-S PFT have to unfold their thin, short-lived leaves at the start of 1082 | THONICKE et al. the growing season, have a disadvantage under this climate. derived from Earth observation data, and allow large-scale estimates Nevertheless over a wide range within the boreal zone, B-NL of FD across Pan-Europe. We demonstrate how FR, FDv and FE vary trees are accompanied by BL-S trees. This means, in this range strongly across Pan-Europe and emerge from environmental and environmental and competitive filtering are not strong enough competitive filtering alike. Co-existence of functionally diverse trees to select for B-NL trees only and both growing strategies can results from plant trait diversity, tree demography and disturbance co-exist. Competitive filtering though pushes BL-S towards lower under varying strength of environmental and competitive filtering. LL and higher WD as in the temperate zone, which minimizes the overlap of trait values between the two PFTs (Figure 1). AC K N OW L E D G M E N T S Due to the BL-S and B-NL co-existence, a wide trait space The authors thank three anonymous reviewers for their careful and (high FR) is covered in which trait combinations are mostly constructive comments that greatly improved the manuscript. This evenly spaced (high FE) and divergent (high FDv), allowing for project has received funding from the European Union's Horizon relatively high resource-use efficiency. The same effect applies 2020 Research and Innovation Programme under grant agreement No further South towards the temperate zone, where B-NL trees 641762 (ECOPOTENTIAL). The study has been supported by the TRY can still be found in mixed temperate forests. initiative on plant traits (http://www.try-db.org). The TRY initiative 2. In the temperate zone, climate supports a longer growing sea- and database is hosted, developed and maintained by J. Kattge and G. son, higher productivity and biomass. Here, BL-S trees dominate Boenisch (MaxPlanck Institute for Biogeochemistry, Jena, Germany) due to competitive filtering. Their tolerance of warmer tempera- in collaboration with the German Centre for Integrative Biodiversity tures together with low carbon investment in short-lived and thin Research (iDiv). We thank Giovanni Gligora for contributing his leaves (Table 1) makes this PFT very competitive along a large data to TRY. TRY is been supported by DIVERSITAS core project climate gradient. Moreover, towards warmer and drier climate bioDISCOVERY, the IGBP, the Global Land Project, the UK Natural conditions, the WD of BL-S trees increases to survive this climatic Environment Research Council (NERC) through its program QUEST stress keeping this PFT dominant. The clear PFT dominance de- (Quantifying and Understanding the Earth System) and the French creases FR, and leads to intermediate FE and FDv. Towards the programmes 'Climat, Environnement, Societe' and 'Fondation pour la Mediterranean zone new tree growth strategies (BL-E, T-NL) Recherche sur la Biodiversite'. The contribution of M. E. Schaepman enter the forest community with milder winter and drier summer is supported by the University of Zurich Research Priority Program changing environmental filtering effects. on Global Change and Biodiversity (URPP GCB). The research car- 3. In the Mediterranean zone, forests are stressed by drought and ried out at the Jet Propulsion Laboratory, California Institute of warm temperatures, therefore, limitations of growing conditions Technology, was under a contract with the National Aeronautics and vary and environmental filtering plays an important role. BL-E and Space Administration. Government sponsorship is acknowledged. T-NL trees with low SLA and high WD are adapted to season- ally dry summers and can therefore co-exist in complementary DATA AVA I L A B I L I T Y S TAT E M E N T and regularly spaced niches which confirms previous findings We have permanently uploaded all relevant data to the public data (Carnicer, Barbeta, Sperlich, Coll, & Penuelas, 2013). Whereas server at the German Research Centre for Geosciences (GFZ) Potsdam, BL-E trees tolerate warm temperatures and drought, T-NL trees Germany, to make the research fully reproducible (Billing et al., 2019). are specialists for warm, dry and fire-prone environments (Table Supporting information link: https​://doi.org/10.5880/PIK.2019.022 1). Given those dynamic changes in niche occupation, FE and FDv are high in mountainous areas and in the semi-arid southern ORCID Mediterranean areas. Where T-NL dominate FDv is slightly lower, Kirsten Thonicke https://orcid.org/0000-0001-5283-4937 that is, niche occupation being less diverse. Werner von Bloh https://orcid.org/0000-0002-7399-2704 Ülo Niinemets https://orcid.org/0000-0002-3078-2192 Fabian D. Schneider https://orcid.org/0000-0003-1791-2041 5 | CO N C LU S I O N REFERENCES Benavides, R., Scherer-Lorenzen, M., & Valladares, F. (2019). The func- We introduce a new generation, large-scale vegetation model for tional trait space of tree species is influenced by the species richness of the canopy and the type of forest. Oikos, 128, 1435–1445. https​:// Europe allowing to dynamically simulate FD on a continental scale, doi.org/10.1111/oik.06348​ alongside with productivity and tree demography. Approximating Bernard-Verdier, M., Navas, M.-L., Vellend, M., Violle, C., Fayolle, A., competition of individual trees with randomly selected functional Garnier, E., & Cornelissen, H. (2012). Community assembly along a trait combinations has proven successful to reproduce trait pat- soil depth gradient: Contrasting patterns of plant trait convergence and divergence in a Mediterranean rangeland. Journal of Ecology, 100, terns, productivity and tree demography of natural forest ecosys- 1422–1433. https​://doi.org/10.1111/1365-2745.12003​ tems. Whereas some trait ranges are constrained by PFTs (SLA), Billing, M., von Bloh, W., Sakschewski, B., & Thonicke, K. (2019) LPJmL- others follow similar climatic gradients across all PFTs (WD). These FIT in Europe. vol 1.0. GFZ Data Services https​://doi.org/10.5880/ results complement the currently available data of trait distributions PIK.2019.022, GFZ Potsdam. THONICKE et al. | 1083 Bohn, U., Zazanashvili, N., & Nakhutsrishvili, G. (2007). The map of the constraints on land surface phenology. Environmental Research natural vegetation of Europe and its application in the Caucasus Letters, 13, 024025. https​://doi.org/10.1088/1748-9326/aaa17b Ecoregion. Bulletin of the Georgia Academy of Science, 175, 112–121. Geller, G. N., Halpin, P. N., Helmuth, B., Hestir, E. L., Skidmore, A., Bruelheide, H., Dengler, J., Purschke, O., Lenoir, J., Jiménez-Alfaro, B., Abrams, M. J., … Williams, K. (2017). Remote sensing for biodiversity. Hennekens, S. M., … Jandt, U. (2018). Global trait-environment rela- In M. Walters, & R. J. Scholes (Eds.), The GEO handbook on biodiver- tionships of plant communities. Nature Ecology & Evolution, 2, 1906– sity observation networks (pp. 187–210). Cham, Switzerland: Springer 1917. https​://doi.org/10.1038/s41559-018-0699-8 Nature. Butler, E. E., Datta, A., Flores-Moreno, H., Chen, M., Wythers, K. R., Hallik, L., Niinemets, Ü., & Wright, I. (2009). Are species shade and drought Fazayeli, F., … Reich, P. B. (2017). Mapping local and global variabil- tolerance reflected in leaf-level structural and functional differentiation ity in plant trait distributions. Proceedings of the National academy of in Northern Hemisphere temperate woody flora? The New Phytologist, Sciences of the United States of America, 114, E10937–E10946. https​:// 184, 257–274. https​://doi.org/10.1111/j.1469-8137.2009.02918.x doi.org/10.1073/pnas.17089​8 4114​ Healey, S. P., Hernandez, M. W., Edwards, D. P., Lefsky, M. A., Freeman, Cadotte, M. W., Carscadden, K., & Mirotchnick, N. (2011). Beyond spe- E., Patterson, P. L., … Lister, A. J. (2015). CMS: GLAS LiDAR-derived cies: Functional diversity and the maintenance of ecological processes Global Estimates of Forest Canopy Height, 2004-2008. Oak Ridge, and services. Journal of Applied Ecology, 48, 1079–1087. https​:// Tennessee, USA. https​://doi.org/10.3334/ORNLD​A AC/1271 doi.org/10.1111/j.1365-2664.2011.02048.x Hisano, M., Searle, E. B., & Chen, H. Y. H. (2018). Biodiversity as a solution Carmona, C. P., de Bello, F., Mason, N. W. H., & Leps, J. (2016). Traits to mitigate climate change impacts on the functioning of forest eco- without borders: Integrating functional diversity across scales. systems. Biological Reviews, 93, 439–456. https​://doi.org/10.1111/ Trends in Ecology & Evolution, 31, 382–394. https​://doi.org/10.1016/j. brv.12351​ tree.2016.02.003 Hooper, D. U., Chapin, F. S., Ewel, J. J., Hector, A., Inchausti, P., Lavorel, Carnicer, J., Barbeta, A., Sperlich, D., Coll, M., & Penuelas, J. (2013). S., … Wardle, D. A. (2005). Effects of biodiversity on ecosystem func- Contrasting trait syndromes in angiosperms and conifers are asso- tioning: A consensus of current knowledge. Ecological Monographs, ciated with different responses of tree growth to temperature on a 75, 3–35. https​://doi.org/10.1890/04-0922 large scale. Frontiers in Plant Science, 4, 409. https​://doi.org/10.3389/ Jetz, W., Cavender-Bares, J., Pavlick, R., Schimel, D., Davis, F. W., Asner, fpls.2013.00409​ G. P., … Ustin, S. L. (2016). Monitoring plant functional diversity Chave, J., Coomes, D., Jansen, S., Lewis, S. L., Swenson, N. G., from space. Nature Plants, 2, 16024. https​://doi.org/10.1038/nplan​ & Zanne, A. E. (2009). Towards a worldwide wood eco- ts.2016.24 nomics spectrum. Ecology Letters, 12, 351–366. https​://doi. Kattge, J., Díaz, S., Lavorel, S., Prentice, I. C., Leadley, P., Bönisch, org/10.1111/j.1461-0248.2009.01285.x G., … Wirth, C. (2011). TRY - A global database of plant Díaz, S., Kattge, J., Cornelissen, J. H. C., Wright, I. J., Lavorel, S., Dray, traits. Global Change Biology, 17, 2905–2935. https​://doi. S., … Gorné, L. D. (2016). The global spectrum of plant form and org/10.1111/j.1365-2486.2011.02451.x function. Nature, 529, 167–171. https​://doi.org/10.1038/natur​ King, D. A., Davies, S. J., Tan, S., & Noor, N. S. M. (2006). The role of e16489 wood density and stem support costs in the growth and mortal- D'Odorico, P., Gonsamo, A., Pinty, B., Gobron, N., Coops, N., Mendez, E., & ity of tropical trees. Journal of Ecology, 94, 670–680. https​://doi. Schaepman, M. E. (2014). Intercomparison of fraction of absorbed pho- org/10.1111/j.1365-2745.2006.01112.x tosynthetically active radiation products derived from satellite data Kunstler, G., Falster, D., Coomes, D. A., Hui, F., Kooyman, R. M., Laughlin, over Europe. Remote Sensing of Environment, 142, 141–154. https​:// D. C., … Westoby, M. (2016). Plant functional traits have globally doi.org/10.1016/j.rse.2013.12.005 consistent effects on competition. Nature, 529, 204–207. https​://doi. Douma, J. C., Witte, J.-P., Aerts, R., Bartholomeus, R. P., Ordoñez, J. C., org/10.1038/natur​e16476 Venterink, H. O., … van Bodegom, P. M. (2012). Towards a functional Laanisto, L., & Niinemets, Ü. (2015). Polytolerance to abiotic stresses: basis for predicting vegetation patterns; incorporating plant traits How universal is the shade-drought tolerance trade-off in woody in habitat distribution models. Ecography, 35, 294–305. https​://doi. species? Global Ecology & Biogeography, 24, 571–580. https​://doi. org/10.1111/j.1600-0587.2011.07140.x org/10.1111/geb.12288​ Ellis, E. C., Kaplan, J. O., Fuller, D. Q., Vavrus, S., Goldewijk, K. K., & Laliberté, E., Legendre, P., & Shipley, B. (2014) FD: Measuring functional Verburg, P. H. (2013). Used planet: A global history. Proceedings of diversity (FD) from multiple traits, and other tools for functional ecology. the National Academy of Sciences of the United States of America, 110, R package version 1.0-12. 7978–7985. https​://doi.org/10.1073/pnas.12172​41110​ Langan, L., Higgins, S. I., & Scheiter, S. (2017). Climate-biomes, pedo-bi- FAO. (2018). Terms and definitions FRA 2020. In D. P. B. A. Pekkarinen omes or pyro-biomes: Which world view explains the tropical for- (Ed.), Forest resources assessment working papers. Rome, Italy: FAO est–savanna boundary in South America? Journal of Biogeography, 44, Forestry Department. 2319–2330. https​://doi.org/10.1111/jbi.13018​ Forkel, M., Carvalhais, N., Schaphoff, S., von Bloh, W., Migliavacca, M., Lausch, A., Bastian, O., Klotz, S., Leitão, P. J., Jung, A., Rocchini, D., … Thurner, M., & Thonicke, K. (2014). Identifying environmental con- Knapp, S. (2018). Understanding and assessing vegetation health by in trols on vegetation greenness phenology through model-data inte- situ species and remote-sensing approaches. Methods in Ecology and gration. Biogeosciences, 11, 7025–7050. https​://doi.org/10.5194/ Evolution, 9, 1799–1809. https​://doi.org/10.1111/2041-210X.13025​ bg-11-7025-2014 Lavorel, S., & Garnier, E. (2002a). Predicting changes in community Funk, J. L., Larson, J. E., Ames, G. M., Butterfield, B. J., Cavender-Bares, composition and ecosystem functioning from plant traits: Revisiting J., Firn, J., … Wright, J. (2017). Revisiting the Holy Grail: Using plant the Holy Grail. Functional Ecology, 16, 545–556. https​://doi. functional traits to understand ecological processes. Biological org/10.1046/j.1365-2435.2002.00664.x Reviews of the Cambridge Philosophical Society, 92, 1156–1173. https​:// Lavorel, S., Diaz, S., Cornelissen, J. H. C., Garnier, E., Harrison, S. P., doi.org/10.1111/brv.12275​ McIntyre, S., … Urcelay, C. (2007). Plant functional types: Are we Garnier, E., Navas, M.-L., & Grigulis, K. (2016). Plant functional diversity. getting any closer to the Holy Grail? In J. G. Canadell, D. E. Pataki, Organisms, traits, community structure, and ecosystem properties. & L. F. Pitelka (Eds.), Terrestrial Ecosystems in a Changing World (pp. Oxford, UK: Oxford University Press. 149–160). Berlin, Heidelberg, New York: Springer. Garonna, I., de Jong, R., Stöckli, R., Schmid, B., Schenkel, D., David, S., Lefsky, M. A., Keller, M., Pang, Y., De Camargo, P. B., & Hunter, M. O. & Schaepman, M. E. (2018). Shifting relative importance of climatic (2007). Revised method for forest canopy height estimation from 1084 | THONICKE et al. Geoscience Laser Altimeter System waveforms. Journal of Applied Prentice, I., Bondeau, A., Cramer, W., Harrison, S., Hickler, T., Lucht, Remote Sensing, 1, 013537-13518. https​://doi.org/10.1117/1.2795724 W., Sykes, M. (2007) Dynamic global vegetation modelling: quan- Liebergesell, M., Reu, B., Stahl, U., Freiberg, M., Welk, E., Kattge, J., … tifying terrestrial ecosystem responses to large-scale environ- Wirth, C. (2016). Functional resilience against climate-driven ex- mental change. In J. G. Canadell, D. E. Pataki, & L. F. Pitelka (Eds.). tinctions - Comparing the functional diversity of European and Terrestrial ecosystems in a changing world. Berlin, Heidelberg, New North American tree floras. PLoS ONE, 11, e0148607. https​://doi. York: Springer. org/10.1371/journ​al.pone.0148607 R-Core-Team. (2019) R: A language and environment for statistical comput- Ma, X., Mahecha, M. D., Migliavacca, M.,van der Plas, F., Benavides, R., ing. R Foundation for Statistical Computing. Retrieved from https​:// Ratcliffe, S., … Wirth, C. (2019). Inferring plant functional diversity www.R-proje​c t.org/. from space: The potential of Sentinel-2. Remote Sensing of Environment, Ratcliffe, S., Liebergesell, M., Ruiz-Benito, P., Madrigal González, J., 233, 111368. https​://doi.org/10.1016/j.rse.2019.111368 Muñoz Castañeda, J. M., Kändler, G., … Wirth, C. (2016). Modes Madrigal-González, J., Ruiz-Benito, P., Ratcliffe, S., Calatayud, J., Kändler, of functional biodiversity control on tree productivity across the G., Lehtonen, A., … Zavala, M. A. (2016). Complementarity effects on European continent. Global Ecology and Biogeography, 25, 251–262. tree growth are contingent on tree size and climatic conditions across https​://doi.org/10.1111/geb.12406​ Europe. Scientific Reports, 6, 10. https​://doi.org/10.1038/srep3​2233 Ratcliffe, S., Wirth, C., Jucker, T.,van der Plas, F., Scherer-Lorenzen, M., Mason, N. W. H., & de Bello, F. (2013). Functional diversity: A tool for Verheyen, K., … Baeten, L. (2017). Biodiversity and ecosystem func- answering challenging ecological questions. Journal of Vegetation tioning relations in European forests depend on environmental Science, 24, 777–780. https​://doi.org/10.1111/jvs.12097​ context. Ecology Letters, 20, 1414–1426. https​://doi.org/10.1111/ Mason, N. W. H., Mouillot, D., Lee, W. G., & Wilson, J. B. (2005). ele.12849​ Functional richness, functional evenness and functional divergence: Reich, P. B. (2014). The world-wide ‘fast–slow’ plant economics spec- The primary components of functional diversity. Oikos, 111, 112– trum: A traits manifesto. Journal of Ecology, 102, 275–301. https​:// 118. https​://doi.org/10.1111/j.0030-1299.2005.13886.x doi.org/10.1111/1365-2745.12211​ Mason, N. W. H., de Bello, F., Mouillot, D., Pavoine, S., & Dray, S. (2013). Reichstein, M., Falge, E., Baldocchi, D., Papale, D., Aubinet, M., Berbigier, A guide for using functional diversity indices to reveal changes in P., … Valentini, R. (2005). On the separation of net ecosystem ex- assembly processes along ecological gradients. Journal of Vegetation change into assimilation and ecosystem respiration: Review and im- Science, 24, 794–806. https​://doi.org/10.1111/jvs.12013​ proved algorithm. Global Change Biology, 11, 1424–1439. https​://doi. MATLAB. (2019). version 9.6 (R2019a). The MathWorks Inc. org/10.1111/j.1365-2486.2005.001002.x Mori, A. S., Lertzman, K. P., & Gustafsson, L. (2017). Biodiversity and Ruiz-Benito, P., Ratcliffe, S., Jump, A. S., Gómez-Aparicio, L., Madrigal- ecosystem services in forest ecosystems: A research agenda for ap- González, J., Wirth, C., … Zavala, M. A. (2017). Functional diver- plied forest ecology. Journal of Applied Ecology, 54, 12–27. https​://doi. sity underlies demographic responses to environmental variation org/10.1111/1365-2664.12669​ in European forests. Global Ecology and Biogeography, 26, 128–141. Mouillot, D., Graham, N. A. J., Villeger, S., Mason, N. W. H., & Bellwood, https​://doi.org/10.1111/geb.12515​ D. R. (2013). A functional approach reveals community responses to Running, S. (2015). MODIS/TERRA Gross Primary Productivity 8-Day L4 disturbances. Trends in Ecology & Evolution, 28, 167–177. https​://doi. Global 500 m SIN Grid V006. Sioux Falls, South Dakota. Retrieved org/10.1016/j.tree.2012.10.004 from https​://lpdaac.usgs.gov. Nachtergaele, F., van Velthuizen, H., Verelst, L., Batjes, N., Dijkshoorn, Sakschewski, B., von Bloh, W., Boit, A., Rammig, A., Kattge, J., Poorter, L., K., van Engelen, V., … Petri, M. (2009). Harmonized world soil database, … Thonicke, K. (2015). Leaf and stem economics spectra drive diver- Food and Agriculture Organization of the United Nations. Retrieved sity of functional plant traits in a dynamic global vegetation model. from http://www.fao.org/soils-porta​l/soil-surve​y/soil-maps-and- Global Change Biology, 21, 2711–2725. https​://doi.org/10.1111/ datab​ases/harmo​nized-world-soil-datab​ase-v12/en/, Rome. gcb.12870​ Naeem, S., Duffy, J. E., & Zavaleta, E. S. (2012). The functions of biological Sandel, B., Gutierrez, A. G., Reich, P. B., Schrodt, F., Dickie, J., & Kattge, diversity in an Age of Extinction. Science, 336, 1401–1406. https​:// J. (2015). Estimating the missing species bias in plant trait mea- doi.org/10.1126/scien​ce.1215855 surements. Journal of Vegetation Science, 26, 828–838. https​://doi. Naveh, Z. (1990) Fire in the Mediterranean - A landscape ecological per- org/10.1111/jvs.12292​ spective. In J. G. Goldammer, & M. J. Jenkins (Eds.), Fire in ecosystem Schaphoff, S., Forkel, M., Müller, C., Knauer, J., von Bloh, W., Gerten, dynamics: Mediterranean and Northern perspectives. The Hague, the D., … Waha, K. (2018). LPJmL4 – A dynamic global vegetation Netherlands: SPB Academic Publ. model with managed land – Part 2: Model evaluation. Geoscientific Niinemets, Ü., & Valladares, F. (2006). Tolerance to shade, drought and Model Development, 11, 1377–1403. https​://doi.org/10.5194/ waterlogging in the temperate dendroflora of the Northern hemi- gmd-11-1377-2018 sphere: Tradeoffs, phylogenetic signal and implications for niche dif- Schaphoff, S., von Bloh, W., Rammig, A., Thonicke, K., Biemans, H., ferentiation. Ecological Monographs, 76, 521–547. Forkel, M., … Waha, K. (2018). LPJmL4 – A dynamic global vegetation Pakeman, R. J. (2011). Functional diversity indices reveal the impacts of land model with managed land – Part 1: Model description. Geoscientific use intensification on plantcommunity assembly. Journal of Ecology, Model Development, 11, 1343–1375. https​://doi.org/10.5194/ 99, 1143–1151. https​://doi.org/10.1111/j.1365-2745.2011.01853.x gmd-11-1343-2018 Papale, D., Reichstein, M., Aubinet, M., Canfora, E., Bernhofer, C., Schneider, F. D., Morsdorf, F., Schmid, B., Petchey, O. L., Hueni, A., Kutsch, W., … Yakir, D. (2006). Towards a standardized processing of Schimel, D. S., & Schaepman, M. E. (2017). Mapping functional di- Net Ecosystem Exchange measured with eddy covariance technique: versity from remotely sensed morphological and physiological for- Algorithms and uncertainty estimation. Biogeosciences, 3, 571–583. est traits. Nature Communications, 8, 1441. https​://doi.org/10.1038/ https​://doi.org/10.5194/bg-3-571-2006 s41467-017-01530-3 Paul-Limoges, E., Damm, A., Hueni, A., Liebisch, F., Eugster, W., Schneider, U., Becker, A., Finger, P., Meyer-Christoffer, A., Rudolf, B., Schaepman, M. E., & Buchmann, N. (2018). Effect of environmen- & Ziese, M. (2011) GPCC Full Data Reanalysis Version 6.0 at 0.5°: tal conditions on sun-induced fluorescence in a mixed forest and a Monthly Land-Surface Precipitation from Rain-Gauges built on GTS- cropland. Remote Sensing of Environment, 219, 310–323. https​://doi. based and Historic Data. https​://doi.org/10.5676/DWD_GPCC/ org/10.1016/j.rse.2018.10.018 FD_M_V6_050. THONICKE et al. | 1085 Šímová, I., Violle, C., Svenning, J.-C., Kattge, J., Engemann, K., Sandel, B., … Enquist, B. J. (2018). Spatial patterns and climate relationships of B I O S K E TC H E S major plant traits in the New World differ between woody and her- baceous species. Journal of Biogeography, 45, 895–916. https​://doi. Kirsten Thonicke is interested in plant trait ecology and its link to org/10.1111/jbi.13171​ ecosystem function and ecosystem services. She approaches this Sitch, S., Smith, B., Prentice, I. C., Arneth, A., Bondeau, A., Cramer, W., by developing and applying a flexible-trait DGVM. She is also in- … Venevsky, S. (2003). Evaluation of ecosystem dynamics, plant ge- terested in developing the concept of functional diversity based ography and terrestrial carbon cycling in the LPJ Dynamic Global on simulated plant traits from the local to the regional scale. She Vegetation Model. Global Change Biology, 9, 161–185. https​://doi. org/10.1046/j.1365-2486.2003.00569.x works on the stabilizing role of biodiversity for ecosystem func- Sitch, S., Huntingford, C., Gedney, N., Levy, P. E., Lomas, M., tioning under climate change. Piao, S. L., … Woodward, F. I. (2008). Evaluation of the terres- trial carbon cycle, future plant geography and climate-carbon Maik Billing works on disentangling the relations among functional cycle feedbacks using five Dynamic Global Vegetation Models (DGVMs). Global Change Biology, 14, 2015–2039. https​://doi. diversity, plant traits, ecosystem functions and their effect on org/10.1111/j.1365-2486.2008.01626.x ecosystem stability. He therefore combines process-based mod- Thonicke, K., Venevsky, S., Sitch, S., & Cramer, W. (2001). The role of fire els and machine learning techniques to further understand what disturbance for global vegetation dynamics: Coupling fire into a dy- makes forests resistant towards the approaching climate crisis. namic global vegetation model. Global Ecology and Biogeography, 10, 661–677. https​://doi.org/10.1046/j.1466-822X.2001.00175.x Thurner, M., Beer, C., Santoro, M., Carvalhais, N., Wutzler, T., Ariane Walz works in the broad field of landscape ecology, on Schepaschenko, D., … Schmullius, C. (2014). Carbon stock and den- impact of climate change to ecosystems, their services and their sity of northern boreal and temperate forests. Global Ecology and adaptive potential. She approaches these questions by ecosys- Biogeography, 23, 297–310. https​://doi.org/10.1111/geb.12125​ Tinner, W., Hubschmid, P., Wehrli, M., Ammann, B., & Conedera, M. tem modelling and seeking solutions that are ecologically as well (1999). Long-term forest fire ecology and dynamics in south- as socially feasible. ern Switzerland. Journal of Ecology, 87, 273–289. https​://doi. org/10.1046/j.1365-2745.1999.00346.x Boris Sakschewski is interested in the adaptive capacity of for- Van Bodegom, P. M., Douma, J. C., Witte, J. P. M., Ordoñez, J. C., Bartholomeus, R. P., & Aerts, R. (2012). Going beyond limita- est ecosystems to environmental change. He approaches this by tions of plant functional types when predicting global ecosys- developing a flexible-trait DGVM which reflects the diversity of tem-atmosphere fluxes: Exploring the merits of traits-based ap- functional plant traits and their plasticity. proaches. Global Ecology and Biogeography, 21, 625–636. https​://doi. org/10.1111/j.1466-8238.2011.00717.x Michael E. Schaepman specializes in Earth Observation using im- Villeger, S., Mason, N. W. H., & Mouillot, D. (2008). New multidi- mensional functional diversity indices for a multifaceted frame- aging spectroscopy combined with forward and inverse physical work in functional ecology. Ecology, 89, 2290–2301. https​://doi. modelling. His interest is to retrieve morphological and physio- org/10.1890/07-1206.1 logical traits from remotely sensed data for biodiversity mapping. Weedon, G. P., Balsamo, G., Bellouin, N., Gomes, S., Best, M. J., & Viterbo, P. (2014). The WFDEI meteorological forcing data set: WATCH Forcing Data methodology applied to ERA-Interim reanal- Author contributions: KTh and MB jointly designed the study and ysis data. Water Resources Research, 50, 7505–7514. https​://doi. jointly conducted the main analyses. KTh, MB, BS, WvB and AW org/10.1002/2014W​R015638 adapted the LPJmL-FIT model to European natural forests. MB Weedon, G. P., Gomes, S., Viterbo, P., Shuttleworth, W. J., Blyth, and WvB performed the simulation experiments. KTh, MB, BS, E., Österle, H., … Best, M. (2011). Creation of the WATCH forc- ing data and its use to assess global and regional reference crop FS and MS analysed the simulated trait data and worked on the evaporation over land during the twentieth century. Journal of computation of the functional diversity indices. KTh wrote the Hydrometeorology, 12, 823–848. https​://doi.org/10.1175/2011J​ first paper draft. ÜN, JP, HC, YO and PvB contributed TRY data to HM1369.1 the study. All authors contributed to the analyses and interpreta- Wright, I. J., Dong, N., Maire, V., Prentice, I. C., Westoby, M., Díaz, S., … Wilf, P. (2017). Global climatic drivers of leaf size. Science, 357, 917– tion of the data, and the final paper draft. 921. https​://doi.org/10.1126/scien​ce.aal4760 Wright, I. J., Reich, P. B., Westoby, M., Ackerly, D. D., Baruch, Z., S U P P O R T I N G I N FO R M AT I O N Bongers, F., … Villar, R. (2004). The worldwide leaf economics Additional supporting information may be found online in the spectrum. Nature, 428, 821–827. https​://doi.org/10.1038/natur​ e02403 Supporting Information section.   Yachi, S., & Loreau, M. (1999). Biodiversity and ecosystem productivity in a fluctuating environment: The insurance hypothesis. PNAS, 96, 1463–1468. How to cite this article: Thonicke K, Billing M, von Bloh W, et Zhang, T., Niinemets, U., Sheffield, J., & Lichstein, J. W. (2018). Shifts al. Simulating functional diversity of European natural forests in tree functional composition amplify the response of forest bio- along climatic gradients. J Biogeogr. 2020;47:1069–1085. mass to climate. Nature, 556, 99–102. https​://doi.org/10.1038/natur​ e26152 https​://doi.org/10.1111/jbi.13809​