Global trait–environment relationships of plant communities

1 Nature Ecology & Evolution December 2018, Volume 2 Issue 12 Pages 1906-1917 https://doi.org/10.1038/s41559-018-0699-8 https://archimer.ifremer.fr/doc/00469/58091/ Archimer https://archimer.ifremer.fr Global trait–environment relationships of plant communities 1, 2, * 2, 3, 4 1, 2 5 Bruelheide Helge , Dengler Jürgen , Purschke Oliver , Lenoir Jonathan , 1, 2, 6 7 8 9 Jiménez-Alfaro Borja , Hennekens Stephan M. , Botta-Dukát Zoltán , Chytrý Milan , 10 11 2, 12 13 10, 12 Field Richard , Jansen Florian , Kattge Jens , Pillar Valério D. , Schrodt Franziska , 2, 12 14 15 16 17 Mahecha Miguel D. , Peet Robert K. , Sandel Brody , Van Bodegom Peter , Altman Jan , 18 19, 20, 21 22 23 Alvarez-Dávila Esteban , Arfin Khan Mohammed A. S. , Attorre Fabio , Aubin Isabelle , 24 25 26 27 Baraloto Christopher , Barroso Jorcely G. , Bauters Marijn , Bergmeier Erwin , Biurrun Idoia, 28 29 30, 31 32 33 Bjorkman Anne D. , Blonder Benjamin , Čarni Andraž , Cayuela Luis , Černý Tomáš , 34 2, 35 36 37 Cornelissen J. Hans C. , Craven Dylan , Dainese Matteo , Derroire Géraldine , 22, 38 39, 40 17, 41 42 De Sanctis Michele , Díaz Sandra , Doležal Jiří , Farfan-Rios William , 43 Feldpausch Ted R. , Fenton Nicole J., Garnier Eric, Guerin Greg R., Gutiérrez Alvaro G., 1, 2 Haider Sylvia , Hattab Tarek, Henry Greg, Hérault Bruno, Higuchi Pedro, Hölzel Norbert, 21 Homeier Jürgen, Jentsch Anke , Jürgens Norbert, Kącki Zygmunt, Karger Dirk N., Kessler Michael, 9 1, 2, 35 Kleyer Michael, Knollová Ilona , Korolyuk Andrey Y., Kühn Ingolf , Laughlin Daniel C., Lens Frederic, Loos Jacqueline, Louault Frédérique, Lyubenova Mariyana I., Malhi Yadvinder, 27 Marcenò Corrado , Mencuccini Maurizio, Müller Jonas V., Munzinger Jérôme, Myers-Smith Isla H., 7 Neill David A., Niinemets Ülo, Orwin Kate H., Ozinga Wim A. , Penuelas Josep, Pérez-Haase Aaron, 17 2 Petřík Petr , Phillips Oliver L., Pärtel Meelis, Reich Peter B., Römermann Christine , 1, 2 1 Rodrigues Arthur V., Sabatini Francesco Maria , Sardans Jordi, Schmidt Marco, Seidler Gunnar , 2 Silva Espejo Javier Eduardo, Silveira Marcos, Smyth Anita, Sporbert Maria , 28 Svenning Jens-Christian , Tang Zhiyao, Thomas Raquel, Tsiripidis Ioannis, Vassilev Kiril, 43 2 1, 2 2 2 Violle Cyrille , Virtanen Risto , Weiher Evan, Welk Erik , Wesche Karsten , Winter Marten , 2, 12 1, 2 Wirth Christian , Jandt Ute 1 Institute of Biology/Geobotany and Botanical Garden, Martin Luther University Halle-Wittenberg, Halle, Germany 2 German Centre for Integrative Biodiversity Research (iDiv), Halle-Jena-Leipzig, Leipzig, Germany 3 Research Group Vegetation Ecology, Institute of Natural Resource Sciences, Zurich University of Applied Sciences, Wädenswil, Switzerland 4 Plant Ecology, Bayreuth Center of Ecology and Environmental Research, University of Bayreuth, Bayreuth, Germany 5 UR ‘Ecologie et Dynamique des Systèmes Anthropisés’ (EDYSAN, UMR 7058 CNRS-UPJV), CNRS, Université de Picardie Jules Verne , Amiens, France 6 Research Unit of Biodiversity (CSIC/UO/PA), University of Oviedo, Campus de Mieres, Mieres, Spain 7 Team Vegetation, Forest and Landscape Ecology, Wageningen Environmental Research (Alterra) , Wageningen, the Netherlands 8 GINOP Sustainable Ecosystems Group, MTA Centre for Ecological Research, Tihany, Hungary 9 Department of Botany and Zoology, Masaryk University, Brno, Czech Republic 10 School of Geography, University of Nottingham, University Park, Nottingham, UK 11 Faculty for Agricultural and Environmental Science, University of Rostock , Rostock, Germany 12 Max Planck Institute for Biogeochemistry, Jena, Germany 13 Department of Ecology, Universidade Federal do Rio Grande do Sul , Porto Alegre, Brazil 14 Department of Biology, University of North Carolina at Chapel Hill , Chapel Hill, NC, USA 15 Department of Biology, Santa Clara University, Santa Clara, CA, USA 16 Department of Conservation Biology, Institute of Environmental Sciences, Leiden University, Leiden, the Netherlands Please note that this is an author-produced PDF of an article accepted for publication following peer review. The definitive publisher-authenticated version is available on the publisher Web site. 2 17 Institute of Botany of the Czech Academy of Sciences, Průhonice, Czech Republic Escuela de Ciencias Agropecuarias y Ambientales – ECAPMA, Universidad Nacional Abierta y a Distancia – UNAD, Sede José Celestino Mutis, Bogotá, Colombia 19 Department of Forestry and Environmental Science, Shahjalal University of Science and Technology , Sylhet, Bangladesh 20 Department of Disturbance Ecology, Bayreuth Center of Ecology and Environmental Research, University of Bayreuth , Bayreuth, Germany 21 Department of Environmental Biology, Sapienza University of Rome , Rome, Italy 22 Great Lakes Forestry Centre, Canadian Forest Service, Natural Resources Canada, Sault Ste Marie, Ontario, Canada 23 Department of Biological Sciences, International Center for Tropical Botany, Florida International University, Miami, FL, USA 24 Campus de Cruzeiro do Su, Universidade Federal do Acre, Acre, Brazil 25 Department of Green Chemistry and Technology (ISOFYS) and Department of Environment (CAVELab), Faculty of Bioscience Engineering, Ghent University, Gent, Belgium 26 Vegetation Analysis & Plant Diversity, Albrecht von Haller Institute of Plant Sciences, University of Göttingen, Göttingen, Germany 27 University of the Basque Country UPV/EHU, Bilbao, Spain 28 Biodiversity Dynamics in a Changing World (BIOCHANGE) & Section for Ecoinformatics & Biodiversity, Department of Bioscience, Aarhus University, Aarhus, Denmark 29 Environmental Change Institute, School of Geography and the Environment, University of Oxford, Oxford, UK 30 Rocky Mountain Biological Laboratory, Crested Butte, CO, USA 31 Institute of Biology, Scientific Research Center of the Slovenian Academy of Sciences and Arts, Ljubljana, Slovenia 32 University of Nova Gorica, Nova Gorica, Slovenia 33 Department of Biology, Geology, Physics and Inorganic Chemistry, Universidad Rey Juan Carlos, Madrid, Spain 34 Department of Forest Ecology, Faculty of Forestry and Wood Science, Czech University of Life Sciences, Prague, Czech Republic 35 Department of Ecological Science, Faculty of Science, Vrije Universiteit Amsterdam, Amsterdam, the Netherlands 36 Department of Community Ecology, Helmholtz Centre for Environmental Research – UFZ , Halle, Germany 37 Department of Animal Ecology and Tropical Biology, University of Würzburg, Würzburg, Germany 38 Cirad, UMR EcoFoG, Campus Agronomique, Kourou, French Guiana 39 Instituto Multidisciplinario de Biología Vegetal, CONICET and FCEFyN, Universidad Nacional de Córdoba, Córdoba, Argentina 40 Department of Biology, Wake Forest University, Winston Salem, NC, USA 41 Herbario Vargas (CUZ), Universidad Nacional de San Antonio Abad del Cusco, Cusco, Peru 42 Geography, College of Life and Environmental Sciences, University of Exeter, Exeter, UK 43 Institut de recherche sur les forêts, Université du Québec en Abitibi-Témiscamingue, Rouyn-Noranda, Quebec, Canada 44 Centre d’Ecologie Fonctionnelle et Evolutive (UMR5175), CNRS, Université de Montpellier, Université Paul-Valéry Montpellier, EPHE, Montpellier, France 45 Terrestrial Ecosystem Research Network, School of Biological Sciences, University of Adelaide, Adelaide, South Australia, Australia 46 Departamento de Ciencias Ambientales y Recursos Naturales Renovables, Facultad de Ciencias Agronómicas, Universidad de Chile, Santiago, Chile 47 UMR 248 MARBEC (CNRS, IFREMER, IRD, UM), Institut Français de Recherche pour l’Exploitation de la MER, Sète, France 48 The Department of Geography, University of British Columbia, Vancouver, British Columbia, Canada 49 Institut National Polytechnique Félix Houphouët-Boigny, Yamoussoukro, Côte d’Ivoire 50 UR Forests & Societies, Cirad, University of Montpellier, Montpellier, France 51 Departamento de Engenharia Florestal, Universidade do Estado de Santa Catarina, Lages, Brazil 52 Institute of Landscape Ecology, University of Münster, Münster, Germany 18 Please note that this is an author-produced PDF of an article accepted for publication following peer review. The definitive publisher-authenticated version is available on the publisher Web site. 3 53 Plant Ecology and Ecosystems Research, University of Göttingen, Göttingen, Germany Biodiversity, Biocenter Klein Flottbek and Botanical Garden, University of Hamburg, Hamburg, Germany 55 Department of Vegetation Ecology, Institute of Environmental Biology, University of Wroclaw, Wrocław, Poland 56 Department of Systematic and Evolutionary Botany, University of Zurich, Zurich, Switzerland 57 Swiss Federal Research Institute WSL, Birmensdorf, Switzerland 58 Landscape Ecology Group, Institute of Biology and Environmental Sciences, University of Oldenburg, Oldenburg, Germany 59 Central Siberian Botanical Garden SB RAS, Novosibirsk, Russia 60 School of Science, Environmental Research Institute, University of Waikato, Hamilton, New Zealand 61 Department of Botany, University of Wyoming, Laramie, WY, USA 62 Naturalis Biodiversity Cente, Leiden University, Leiden, the Netherlands 63 Agroecology, University of Göttingen, Göttingen, Germany 64 UCA, INRA, VetAgro Sup, UREP, Clermont-Ferrand, France 65 Department of Ecology and Environmental Protection, Faculty of Biology, University of Sofia, Sofia, Bulgaria 66 Environmental Change Institute, School of Geography and the Environment, University of Oxford, Oxford, UK 67 ICREA, Barcelona, Spain 68 CREAF, Barcelona, Spain 69 Millennium Seed Bank, Conservation Science, Royal Botanic Gardens Kew, Ardingly, UK 70 AMAP, IRD, CIRAD, CNRS, INRA, Université Montpellier, Montpellier, France 71 School of GeoSciences, University of Edinburgh, Edinburgh, UK 72 Conservación y Manejo de Vida Silvestre, Universidad Estatal Amazónica, Puyo, Ecuador 73 Department of Crop Science and Plant Biology, Estonian University of Life Science, Tartu, Estonia 74 Landcare Research, Lincoln, New Zealand 75 Institute for Water and Wetland Research, Radboud University Nijmegen, Nijmegen, the Netherlands 76 Global Ecology Unit, CREAF-CEAB-UAB, CSIC, Cerdanyola del Vallès, Spain 77 Department of Evolutionary Biology, Faculty of Biology, Ecology and Environmental Sciences, University of Barcelona, Barcelona, Spain 78 Spanish Research Council (CEAB-CSIC), Center for Advanced Studies of Blanes, Blanes, Spain 79 School of Geography, University of Leeds, Leeds, UK 80 University of Tartu, Tartu, Estonia 81 Department of Forest Resources, University of Minnesota, St. Paul, MN, USA 82 Hawkesbury Institute for the Environment, Western Sydney University, Sydney, New South Wales, Australia 83 Institute of Ecology and Evolution, Friedrich Schiller University Jena, Jena, Germany 84 Departamento de Engenharia Florestal, Universidade Regional de Blumenau, Blumenau, Brazil 85 Data and Modelling Centre, Senckenberg Biodiversity and Climate Research Centre (BiK-F), Frankfurt am Main, Germany 86 Department of Biology University of La Serena, La Serena, Chile 87 Laboratório de Botânica e Ecologia Vegetal, Centro de Ciências Biológicas e da Natureza, Museu Universitário, Universidade Federal do Acre, Rio Branco, Brazil 88 College of Urban and Environmental Sciences, Peking University, Beijing, China 89 Iwokrama International Centre for Rain Forest Conservation and Development, Georgetown, Guyana 90 Department of Botany, School of Biology, Aristotle University of Thessaloniki, Thessaloniki, Greece 91 Institute of Biodiversity and Ecosystem Research, Bulgarian Academy of Sciences, Sofia, Bulgaria 92 Department of Physiological Diversity, Helmholtz Center for Environmental Research – UFZ, Leipzig, Germany 93 Department of Ecology & Genetics, University of Oulu, Oulu, Finland 94 Department of Biology, University of Wisconsin – Eau Claire, Eau Claire, WI, USA 95 Senckenberg Museum of Natural History Görlitz, Görlitz, Germany 96 International Institute (IHI) Zittau, TU Dresden, Zittau, Germany 97 Systematic Botany and Functional Biodiversity, University of Leipzig, Leipzig, Germany 54 Please note that this is an author-produced PDF of an article accepted for publication following peer review. The definitive publisher-authenticated version is available on the publisher Web site. 4 * Corresponding author : Helge Bruelheide Abstract : Plant functional traits directly affect ecosystem functions. At the species level, trait combinations depend on trade-offs representing different ecological strategies, but at the community level trait combinations are expected to be decoupled from these trade-offs because different strategies can facilitate coexistence within communities. A key question is to what extent community-level trait composition is globally filtered and how well it is related to global versus local environmental drivers. Here, we perform a global, plot-level analysis of trait–environment relationships, using a database with more than 1.1 million vegetation plots and 26,632 plant species with trait information. Although we found a strong filtering of 17 functional traits, similar climate and soil conditions support communities differing greatly in mean trait values. The two main community trait axes that capture half of the global trait variation (plant stature and resource acquisitiveness) reflect the trade-offs at the species level but are weakly associated with climate and soil conditions at the global scale. Similarly, within-plot trait variation does not vary systematically with macro-environment. Our results indicate that, at fine spatial grain, macroenvironmental drivers are much less important for functional trait composition than has been assumed from floristic analyses restricted to co-occurrence in large grid cells. Instead, trait combinations seem to be predominantly filtered by local-scale factors such as disturbance, fine-scale soil conditions, niche partitioning and biotic interactions. Please note that this is an author-produced PDF of an article accepted for publication following peer review. The definitive publisher-authenticated version is available on the publisher Web site. 233 Introduction 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 How climate drives the functional characteristics of vegetation across the globe has been a key question in ecological research for more than a century1. While functional information is available for a large portion of the global pool of plant species, we do not know how functional traits of the different species that co-occur in a community are combined, which is what determines their joint effect on ecosystems2-4. At the species level, Díaz et al.5 demonstrated that 74% of the global spectrum of six key plant traits determining plant fitness in terms of survival, growth and reproduction can be accounted for by two principal components (PCs). They showed that the functional space occupied by vascular plant species is strongly constrained by trade-offs between traits and converges on a small set of successful trait combinations, confirming previous findings6-9. While these constraints describe evolutionarily viable ecological strategies for vascular plant species globally, they provide only limited insight into trait composition within communities. There are many reasons why trait composition within communities would produce very different patterns, and indeed much theory predicts this10-11. However, it is still unknown to what extent community-level trait composition depends on local factors (microclimate, fine-scale soil properties, disturbance regime10, successional dynamics2) and regional to global environmental drivers (macroclimate6,12-13, coarse-scale soil properties3,14). As ecosystem functions and services are ultimately dependent on the traits of the species composing ecological communities, exploring community trait composition at the global scale can advance our understanding of how climate change and other anthropogenic drivers affect ecosystem functioning. 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 So far, studies relating trait composition to the environment at continental to global extents have been restricted to coarse-grained species occurrence data (e.g. presence in 1° grid cells1517 ). Such data capture neither biotic interactions (co-occurrence in large grid cells does not indicate local co-existence), nor local variation in environmental filters (e.g. variation in soil, topography or disturbance regime within grid cells). In contrast, functional composition of ecological communities sampled at fine-grained vegetation plots – with areas of few to a few hundred square meters – is the direct outcome of the interaction between both local and largescale factors. Here, we present a global analysis of plot-level trait composition. We combined the ‘sPlot’ database, a new global initiative incorporating more than 1.1 million vegetation plots from over 100 databases (mainly forests and grasslands; see Methods), with 30 largescale environmental variables and 18 key plant functional traits derived from TRY, a global plant-trait database (see Methods, Table 2). We selected these 18 traits because they affect different key ecosystem processes and are expected to respond to macroclimatic drivers (Table 1). In addition, they were sufficiently measured across all species globally to allow for imputation of missing values (see Methods). All analyses were confined to vascular plant species and included all vegetation layers in a community, from the canopy to the herb layer (see Methods). 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 We used this unprecedented fine-resolution dataset to test the hypothesis (Hypothesis 1) that plant communities show evidence of environmental or biotic filtering at the global scale, making use of the observed variation of plot-level trait means and means of within-plot trait variation across communities. Ecological theory suggests that community-level convergence could be interpreted as the result of filtering processes, including environmental filtering and biotic interactions. Globally, temperature and precipitation drive the differences in vegetation between biomes, suggesting strong environmental filtering3,11 that constrains the number of successful trait combinations and leads to community-level trait convergence. Similarly, biotic interactions may eliminate excessively divergent trait combinations18,19. However, alternative functional trait combinations may confer equal fitness in the same environment10. If plant communities show a global variation of plot-level trait means higher than expected by chance, and a lower than expected within-plot trait variation (see Figure 1), this would support the view that environmental or biotic filtering are dominant structuring processes of community trait composition at the global scale. A consequence of strong community-level trait convergence, and thus low variation within plots with species trait values centred around the mean, would be that plot-level means will be similar to the trait values of the species in that plot. Hence, community mean trait values should then mirror the trait values of individual species5. 289 290 291 292 293 294 295 296 While Hypothesis 1 addresses the degree of filtering, it does not make a statement on the attribution of driving factors. The main drivers should correlate strongly (though not necessarily linearly20) with plot-level trait means and within-plot trait variance. Identifying these drivers has the potential to fundamentally improve our understanding of global traitenvironment relationships. We tested the hypothesis (Hypothesis 2) that there are strong correlations between global environmental drivers such as macroclimate and coarse-scale soil properties and both plot-level trait means and within-plot trait variances3,6,12-17,20-24 (see Table 1 for expected relationships and Supplementary Table 2 for variables used). Such evidence, 297 298 although correlative, may contribute to the formulation of novel hypotheses to explain global plant trait patterns. 299 300 Results and Discussion 301 302 303 304 305 306 307 308 Consistent with Hypothesis 1 and as illustrated in Figure 1, global variation in plot-level trait means was much higher than expected by chance: all traits had positive standardized effect sizes (SESs), which were significantly > 0 for 17 out of 18 traits based on gap-filled data (mean SES = 8.06 standard deviations (SD), Table 2). This suggests that environmental or biotic filtering is a dominant force of community trait composition globally. Also as predicted by Hypothesis 1, within-plot trait variance was typically lower than expected by chance (mean SES = -1.76 SD, significantly < 0 for ten traits but significantly > 0 for three traits; Table 2). Thus, trait variation within communities may also be constrained by filtering. 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 Trait correlations at the community level were relatively well captured by the first two axes of a Principal Component Analysis (PCA) for both plot-level trait means and within-plot trait variances (Figures 1 and 2). The dominant axes were determined by those traits with the highest absolute SESs of plot-level trait mean trait values (Table 2, mean of CWMs). The PCA of plot-level trait means (Fig. 2) reflects two main functional continua on which community trait values converge: one from short-stature, small-seeded communities such as grasslands or herbaceous vegetation to tall-stature communities with large, heavy diaspores such as forests (the size spectrum), and the other from communities with resource-acquisitive to those with resource-conservative leaves (i.e. the leaf economics spectrum)7. The high similarity between this PCA and the one at the species level by Díaz et al.5 is striking: here at the community level, based on 1.1 million plots, the same functional continua emerged as at the species level, based on 2,214 species. While the trade-offs between different traits at the species level can be understood from a physiological and evolutionary perspective, finding similar trade-offs between traits at the community level was unexpected, as species with opposing trait values can co-exist in the same community. In combination with our finding of strong trait convergence, these results reveal a strong parallel of present-day community assembly to individual species’ evolutionary histories. 326 327 328 329 330 331 332 333 334 335 336 Surprisingly, we found only limited support for Hypothesis 2. Community-level trait composition was poorly captured by global climate and soil variables. None of the 30 environmental variables accounted individually for more than 10% of the variance in the traits defining the main dimensions in Fig. 2 (Supplementary Fig. 2). The coefficients of determination were not improved when testing for non-linear relationships (see Methods). Using all 30 environmental variables simultaneously as predictors only accounted for 10.8% or 14.0% of the overall variation in plot-level trait means (cumulative variance, respectively, of the first two or all 18 constrained axes in a Redundancy Analysis). Overall, our results show that similar global-scale climate and soil conditions can support communities that differ markedly in mean trait values and that different climates can support communities with rather similar mean trait values. 337 338 339 340 341 342 343 344 345 346 347 348 The ordination of within-plot variance of the different traits (Fig. 3) revealed two main continua. Variances of plant height and diaspore mass varied largely independently of variances of traits representing the leaf economics spectrum. This suggests that short and tall species can be assembled together in the same community independently from combining species with acquisitive leaves with species with conservative leaves. Global climate and soil variables accounted for even less variation on the first two PCA axes in within-plot trait variances than on the first two PCA axes in plot-level trait means. Only two environmental variables had r2 > 3% (Supplementary Fig. 3), whether allowing for non-linear relationships (see Methods) or not, and overall, macro-environment accounted for only 3.6% or 5.0% of the variation (cumulative variance, respectively, of the first two or all 18 constrained axes). Removing species richness effects from within-plot trait variances did not increase the amount of variation explained by the environment (see Methods). 349 350 351 352 353 354 355 356 357 358 359 The findings of our study contrast strongly with studies where the variation in traits between species was calculated at the level of the species pool in large grid cells15,16, suggesting that plot-level and grid cell-level trait composition are driven by different factors21. Plot-level trait means and variances may both be predominantly driven by local environmental factors, such as topography (e.g. north- vs. south-facing slopes), local soil characteristics (e.g. soil depth and nutrient supply)3,14,24,25, disturbance regime (including land use26 and successional status2,27) or biotic interactions18-19,28, while broad-scale climate and soil conditions may only become relevant for the whole species pool in large grid cells. Such differences emphasize the importance of local environment in affecting the communities’ trait composition and should be taken into account when interpreting the effect of environmental drivers in functional trait diversity using data on either floristic pools or ecological communities. 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 We note that the strongest community-level correlations with environment were found for traits not linked to the leaf economics spectrum. Mean stem specific density increased with potential evapotranspiration (PET, r2=15.6%; Fig. 4a, b), reflecting the need to produce denser wood with increasing evaporative demand. Leaf N:P ratio increased with growingseason warmth (growing degree days above 5°C, GDD5, r2=11.5%; Fig. 4d), indicating strong phosphorus limitation29 in most plots in the tropics and subtropics (Fig. 4c, d). This pattern was not brought about by a parallel increase in the presence of legumes, which tend to have relatively high N:P ratios; excluding all species of Fabaceae resulted in a very similar relationship with GDD5 (r2=10.0%). The global N:P pattern is consistent with results based on traits of single species related to mean annual temperature30. We assume that the main underlying mechanism is the high soil weathering rate at high temperatures and humidity, which in the tropics and subtropics was not reset by Pleistocene glaciation. Thus, phosphorus limitation may weaken the relationships between productivity-related traits and macroclimate (Supplementary Fig. 2). For example, specific leaf area may be low as consequence of low nutrient availability3,14,24-25 in favourable climates as well as be low as consequence of low temperature and precipitation under favourable nutrient supply. Overall, our findings are relevant in improving Dynamic Global Vegetation Models (DGVMs), which so far have used trait information only from a few calibration plots22. The sPlot database provides muchneeded empirical data on the community trait pool in DGVMs31 and identifies traits that 379 380 should be considered when predicting ecosystem functions from vegetation, such as stem specific density and leaf N:P ratio. 381 382 383 384 385 386 387 388 389 390 391 392 393 394 Our results were surprisingly robust both to the selection of trait data, when comparing different plant formations and when explicitly accounting for the uneven distribution of plots. Using the original trait values measured for the species from the TRY database for the six traits used by Díaz et al.5 (see Methods), resulted in the same two main functional continua and an overall highly similar ordination pattern (Supplementary Fig. 4) compared to using gap-filled data for 18 traits (Fig. 2). Community-level trait composition was also similarly poorly captured by global climate and soil variables. Single regressions of CWMs with all environmental variables revealed very similar patterns to those based on gap-filled traits (Supplementary Fig. 5). Similarly, subjecting the CWMs based on six original traits to a Redundancy Analysis with all 30 environmental variables accounted only for 20.6% or 21.8% of the overall variation in CWMs (cumulative variance of the first two or all six constrained axes, respectively, Supplementary Fig. 4). These results clearly demonstrate that the imputation of missing trait values did not result in spurious artefacts which may have obscured community trait-environment relationships. 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 We also assessed whether the observed trait-environment relationships hold for forests and non-forest vegetation independently (see Methods). Both subsets confirmed the overall patterns in trait means (Supplementary Figs. 3-6). The variance in plot-level trait means explained by large-scale climate and soil variables was higher for forest than non-forest plots, probably because forests belong to a well-defined and rather resource-conservative formation, whereas non-forest plots encompass a heterogeneous mixture of different vegetation types, ranging from alpine meadows to semi-deserts, and tend to depend more on disturbance and management, which can strongly affect trait-environment relationships of communities21. Finally, to test whether our findings depended on the uneven distribution of plots among the world’s different climates and soils, we repeated the analyses in 100 subsets of ~100,000 plots resampled in the global climate space (Supplementary Figs. 7-8). The analyses of the resampled datasets revealed the same patterns and confirmed the impact of PET and GDD5 on stem specific density and leaf N:P ratio, respectively. The correlations between trait means and environmental variables were, however, stronger in the resampled subsets, possibly because the resampling procedure reduced the overrepresentation of the temperate-zone areas with intermediate climatic values. 411 412 413 414 415 416 417 418 419 420 421 Our findings have important implications for understanding and predicting plant community trait assembly. First, worldwide trait variation of plant communities is captured by a few main dimensions of variation, which are surprisingly similar to those reported by species-based studies5,7-9, suggesting that the drivers of past trait evolution, which resulted in the present-day species-level trait spectra5, are also reflected in the composition of today’s plant communities. If species-level trade-offs indeed constrain community assembly, then the present-day contrasts in trait composition of terrestrial plant communities should also have existed in the past and will probably remain, even for novel communities, in the future. Most species in our present-day communities evolved under very variable filtering conditions across the globe, with respect to temperature and precipitation regimes. Therefore, it can be assumed that future filtering conditions will result in novel communities that follow the same functional continua 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 from short-stature, small-seeded communities to tall-stature communities with large, heavy diaspores and from communities with resource-acquisitive to those with resourceconservative leaves. Second, the main plot-level vegetation trait continua cannot easily be captured by coarse-resolution environmental variables21. This brings into question both the use of simple large-scale climate relationships to predict the leaf economics spectra of global vegetation13,15-16,22 and attempts to derive net primary productivity and global carbon and water budgets from global climate, even when employing powerful trait-based vegetation models31. The finding that within-plot trait variances were only very weakly related to global climate or soil variables points to the importance of i) local-scale climate or soil variables, ii) disturbance regimes or iii) biotic interactions for the degree of local trait dispersion11. Finally, both our findings on the limited role of large-scale climate in explaining trait patterns and on the prevalence of phosphorus limitation in most plots in the tropics and subtropics call for including local variables when predicting community trait patterns. Even under similar macro-environmental conditions, communities can vary greatly in trait means and variances, consistent with high local variation in species’ trait values3,6-7. Future research on functional response of communities to changing climate should incorporate the effect of local environmental conditions24-26 and biotic interactions18-19 for building reliable predictions of vegetation dynamics. 440 441 Material and Methods 442 443 444 445 446 447 448 449 Vegetation Data. The sPlot 2.1 vegetation database contains 1,121,244 plots with 23,586,216 species × plot observations, i.e. records of a species in a plot (https://www.idiv.de/en/sdiv/working_groups/wg_pool/splot.html). This database aims at compiling plot-based vegetation data from all vegetation types worldwide, but with a particular focus on forest and grassland vegetation. Although the initial aim of sPlot was to achieve global coverage, the plots are very unevenly distributed with most data coming from Europe, North America and Australia and an overrepresentation of temperate vegetation types (Supplementary Fig. 1). 450 451 452 453 454 455 456 457 458 459 460 461 462 463 For most plots (97.2%) information on the single species’ relative contribution to the sum of plants in the plot was available, expressed as cover, basal area, individual count, importance value or per cent frequency in subplots. For the other 2.8% (31,461 plots), for which only presence/absence (p/a) was available, we assigned equal relative abundance to the species (1/species richness). For plots with a mix of cover and p/a information (mostly forest plots, where herb layer information had been added on a p/a basis; 8,524 plots), relative abundance was calculated by assigning the smallest cover value that occurred in a particular plot to all species with only p/a information in that plot. In most cases (98.4%), plot records in sPlot include full species lists of vascular plants. Bryophytes and lichens were additionally identified in 14% and 7% of plots, respectively. After removing plots without geographic coordinates and all observations on bryophytes and lichens, the database contained 22,195,966 observations on the relative abundance of vascular plant species in a total of 1,117,369 plots. The temporal extent of the data spans from 1885 to 2015, but >95% of vegetation plots were recorded later than 1980. Plot size was reported in 65.4% of plots. 464 465 While forest plots had plot sizes ≥100 m2, and in most cases ≤1,000 m2, non-forest plots typically ranged from 5 to 100 m2. 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 Taxonomy. To standardize the nomenclature of species within and between sPlot and TRY (see below), we constructed a taxonomic backbone of the 121,861 names contained in the two databases. Prior to name matching, we ran a series of string manipulation routines in R, to remove special characters and numbers, as well as standardized abbreviations in names. Taxon names were parsed and resolved using Taxonomic Name Resolution Service version 4.0 (TNRS32; http://tnrs.iplantcollaborative.org; accessed 20 Sep 2015), selecting the best match across the five following sources: i) The Plant List (version 1.1; http://www.theplantlist.org/; Accessed 19 Aug 2015), ii) Global Compositae Checklist (GCC, http://compositae.landcareresearch.co.nz/Default.aspx; accessed 21 Aug 2015), iii) International Legume Database and Information Service (ILDIS, http://www.ildis.org/LegumeWeb; accessed 21 Aug 2015), iv) Tropicos (http://www.tropicos.org/; accessed 19 Dec 2014), and v) USDA Plants Database (http://usda.gov/wps/portal/usda/usdahome; accessed 17 Jan 2015). We allowed for partial matching to the next higher taxonomic rank (genus or family) in cases where full taxon names could not be found. All names matched or converted from a synonym by TNRS were considered accepted taxon names. In cases when no exact match was found (e.g. when alternative spelling corrections were reported), names with probabilities of ≥ 95% or higher were accepted and those with < 95% were examined individually. Remaining non-matching names were resolved based on the National Center for Biotechnology Information's Taxonomy database (NCBI, http://www.ncbi.nlm.nih.gov/; accessed 25 Oct 2011) within TNRS, or sequentially compared directly against The Plant List and Tropicos (accessed September 2015). Names that could not be resolved against any of these lists were left as blanks in the final standardized name field. This resulted in a total of 86,760 resolved names, corresponding to 664 families, occurring in sPlot or TRY or both. Classification into families was carried out according to APGIII33, and was used to identify non-vascular plant species (~5.1% of the taxon names) which were excluded from the subsequent statistical analysis. 492 493 494 495 496 497 498 499 500 501 502 503 504 505 Trait Data. Data for 18 traits that are ecologically relevant (Table 1) and sufficiently covered across species34 were requested from TRY35 (version 3.0) on the 10th August, 2016. We applied gap-filling with Bayesian Hierarchical Probabilistic Matrix Factorization (BHPMF34,36-37). We used the prediction uncertainties provided by BHPMF for each imputation to assess the quality of gap-filling and removed all imputations with a coefficient of variation > 137. We obtained 18 gap-filled traits for 26,632 out of a total of 58,065 taxa in sPlot, which corresponds to 45.9% of all species but to 88.7% of all species × plot combinations. Trait coverage of the most frequent species was 77.2% and 96.2% for taxa that occurred in more than 100 or 1,000 plots, respectively. The gap-filled trait data comprised observed and imputed values on 632,938 individual plants, which we loge transformed and aggregated by taxon. For those taxa that were recorded at the genus level only, we calculated genus means. Out of 22,195,966 records of vascular plant species with geographic reference, 21,172,989 (=95.4%) refer to taxa for which we had gap-filled trait values. This resulted in 1,115,785 and 1,099,463 plots for which we had at least one taxon or two taxa with a trait 506 507 value (99.5% and 98.1%, respectively, of all 1,121,244 plots), and for which trait means and variances could be calculated. 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 As some mean values of traits in TRY were based on a very small number of replicates per species, which results in uncertainty in trait mean and variance calculations38, we tested to which degree the trait patterns in the dataset might be caused by a potential removal of trait variation by imputation of trait values and additionally carried out all analyses using the original trait data on the same 632,938 individual plants instead of gap-filled data (Supplementary Table 1). The degree of trait coverage of species ranged between 7.0% and 58.0% for leaf fresh mass and plant height, respectively. Across all species, mean coverage of species with original trait values was 21.8%, as compared to 45.9% for gap-filled trait data. Linking these trait values to the species occurrence data resulted in a coverage of species × plot observations with trait values between 7.6% and 96.6% for conduit element length and plant height, respectively (Supplementary Table 1), with a mean of 60.7% as compared to 88.7% for those based on gap-filled traits. Using these original trait values to calculate community-weighted mean (CWM) trait values (see below) resulted in a plot coverage of trait values between 48.2% and 100% for conduit element length and SLA, respectively. Across all plots, mean coverage of plots with original trait values was 89.3%, as compared to 100% for gap-filled trait data (Supplementary Table 1). 524 525 526 We are aware that using species mean values for traits excludes the possibility to account for intraspecific variance, which can also strongly respond to the environment39. Thus, using one single value for a species is a source of error in calculating trait means and variances. 527 528 529 530 531 532 533 534 Environmental Data. We compiled 30 environmental variables (Supplementary Table 2). Macroclimate variables were extracted from CHELSA40-41, V1.1 (Climatologies at High Resolution for the Earth’s Land Surface Areas, www.chelsa-climate.org). CHELSA provides 19 bioclimatic variables equivalent to those used in WorldClim (www.worldclim.org) at a resolution of 30 arc sec (~ 1 km at the equator), averaging global climatic data from the period 1979–2013 and using a quasi-mechanistic statistical downscaling of the ERA-Interim reanalysis42. 535 536 537 538 539 540 541 542 Variables reflecting growing-season warmth were growing degree days above 1°C (GDD1) and 5°C (GDD5), calculated from CHELSA data43. We also compiled an index of aridity (AR) and a model for potential evapotranspiration (PET) extracted from the Consortium of Spatial Information (CGIAR-CSI) website (www.cgiar-csi.org). In addition, seven soil variables were extracted from the SOILGRIDS project (https://soilgrids.org/, licensed by ISRIC – World Soil Information), downloaded at 250 m resolution and then resampled using the 30 arc second grid of CHELSA (Supplementary Table 2). We refer to these climate and soil data as “environmental data”. 543 Community trait composition. 544 545 For every trait j and plot k, we calculated the plot-level trait means as community-weighted mean (CWM) according to2,44: , 546 547 548 549 550 , = , = , ( , ( , − , − , ) CWV is equal to functional dispersion as described by Rao’s quadratic entropy46, when using a squared Euclidean distance matrix di,j,k 47: , 553 554 , where nk is the number of species sampled in plot k, pi,k is the relative abundance of species i in plot k, referring to the sum of abundances for all species with traits in the plot, and ti,j is the mean value of species i for trait j. This computation was done for each of the 18 traits for 1,115,785 plots. The within-plot trait variance is given by community-weighted variance (CWV)44,45: , 551 552 = ) = = , , , , We had CWV information for 18 traits for 1,099,463 plots, as at least two taxa were needed to calculate CWV. We performed the calculations using the 'data.table' package48 in R. 555 556 557 558 559 560 561 562 563 Assessing the degree of filtering. To analyse how plot-level trait means and within-plot trait variances (based on gap-filled trait data) depart from random expectation, for each trait we calculated standardized effect sizes (SESs) for the variance in CWMs and for the mean in CWVs. Significantly positive SESs in variance of CWM and significantly negative ones in the mean of CWV can be considered a global-level measure of environmental or biotic filtering. To provide an indication of the global direction of filtering, we also report SESs for the mean of CWM trait values. Similarly, to measure how much within-community trait dispersion varied globally, we also calculated SESs for the variance in CWV. 564 565 566 567 568 569 570 571 572 573 SESs were obtained from 100 runs of randomizing trait values across all species globally. In every run we calculated CWM and CWV with random trait values, but keeping all species abundances in plots. Thus, the results of randomization are independent from species cooccurrences structure of plots49. For every trait, the SESs of the variance in CWM were calculated as the observed value of variance in CWM minus the mean variance in CWM of the random runs, divided by the standard deviation of the variance in CWM of the random runs (Fig. 1). SESs for the mean in CWM, the mean in CWV and the variance in CWV were calculated accordingly. Tests for significance of SESs were obtained by fitting generalized Pareto-distribution of the most extreme random values and then estimating p values form this fitted distribution50. 574 575 576 Vegetation trait-environment relationships. Out of the 1,115,785 plots with CWM values, 1,114,304 (99.9%) had complete environmental information and coordinates. This set of plots 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 was used to calculate single linear regressions of each of the 18 traits on each of the 30 environmental variables. We used the 'corrplot' function51 in R to illustrate Pearson correlation coefficients (see Supplementary Figs. 1-2, 4, 6, 8) and for the strongest relationships produced bivariate graphs and mapped the global distribution of the CWM values using kriging interpolation in ArcGIS 10.2 (Fig. 4). We also tested for non-linear relationships with environment by including an additional quadratic term in the linear model and then report coefficients of determination. As in the linear relationships of CWM with environment, the highest r2 values in models with an additional quadratic term were encountered between stem specific density and PET (r2=0.156) and leaf N:P ratio and growing degree days above 5°C (GDD5, r2=0.118). These were not substantially different from the linear CWM-environment relationships, which had r2=0.156 and r2=0.115, respectively (Fig. 4, Supplementary Fig. 2). Similarly, including a quadratic term in the regressions did not increase the CWV-environment correlations. Here, the strongest correlations were encountered between plant height and soil pH (r2=0.044) and between specific leaf area (SLA) and the volumetric content of coarse fragments in the soil (CoarseFrags, r2=0.037), which were similar to those in the linear regressions (r2=0.029 and r2=0.036, respectively, Supplementary Fig. 3). 594 595 596 597 598 599 600 601 To account for a possible confounding effect of species richness on CWV, which may cause low CWV through competitive exclusion of species, we regressed CWV on species richness and then calculated all Pearson correlation coefficients with the residuals of this relationship against all climatic variables. Here, the highest correlation coefficients were encountered between PET and CWV of conduit element length (r2=0.038), followed by the relationship of specific leaf area (SLA) and the volumetric content of coarse fragments in the soil (CoarseFrags, r2=0.034), which were very similar in magnitude to the CWV environment correlations (r2=0.035 and r2=0.036, respectively; Supplementary Fig. 3). 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 The CWMs and CWVs were scaled to a mean of zero and standard deviation of one and then subjected to a Principal Component Analysis (PCA), calculated with the 'rda' function from the 'vegan' package52. Climate and soil variables were fitted post hoc to the ordination scores of plots of the first two axes, producing correlation vectors using the 'envfit' function. We refrain from presenting any inference statistics, as with > 1.1 million plots all environmental variables showed statistically significant correlations. Instead, we report coefficients of determination (r2), obtained from Redundancy Analysis (RDA), using all 30 environmental variables as constraining matrix, resulting in a maximum of 18 constrained axes corresponding to the 18 traits. We report both r2 values of the first two axes explained by environment, which is the maximum correlation of the best linear combination of environmental variables to explain the CWM or CWV plot × trait matrix and r2 values of all 18 constrained axes explained by environment. We plotted the PCA results using the 'ordiplot' function and coloured the points according to the logarithm of the number of plots that fell into grid cells of 0.002 in PCA units (resulting in approximately 100,000 cells). For further details, see the captions of the figures. 617 618 619 Additionally, we carried out the PCA and RDA analyses, using CWMs based on original trait values (see above). Because of a poor coverage of some traits we confined the analyses with original trait values to the six traits used by Díaz et al.5, which were leaf area, specific leaf 620 621 area, leaf N, seed mass, plant height and stem specific density. Using these six traits resulted in 954,459 plots that had at least one species with a trait value for each of the six traits. 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 Testing for formation-specific patterns. We carried out separate analyses for two ‘formations’: forest and for non-forest plots. We defined as forest plots that had > 25% cover of the tree layer. However, this information was available for only 25% of the plots in our sPlot database. Thus, we also assigned formation status based on growth form data from the TRY database. We defined plots as ‘forest’ if the sum of relative cover of all tree taxa was > 25%, but only if this did not contradict the requirement of > 25% cover of the tree layer (for those records for which this information was given in the header file). Similarly, we defined non-forest plots by calculating the cover of all taxa that were not defined as trees and shrubs (also taken from the TRY plant growth form information) and that were not taller than 2 m, using the TRY data on mean plant height. We assigned the status ‘non-forest’ to all plots that had >90% cover of these low-stature, non-tree and non-shrub taxa. In total, 21,888 taxa out of the 52,032 in TRY which also occurred in sPlot belonged to this category, and 16,244 were classed as trees. The forests and non-forest plots comprised 330,873 (29.7%) and 513,035 (46.0%) of all plots, respectively. We subjected all CWM values for forest and non-forest plots to PCA, RDA and bivariate linear regressions to environmental variables as described above. 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 The forest plots, in particular, confirmed the overall patterns, with respect to variation in CWM explained by the first two PCA axes (60.5%) and the two orthogonal continua from small to large size and the leaf economics spectrum (Supplementary Fig. 6). The variation explained by macroclimate and soil conditions was much larger for the forest subset than for the total data, with the best relationship (leaf N:P ratio and the mean temperature of the coldest quarter, bio11) having r2=0.369 and the second next best ones (leaf N:P ratio and GDD1 and GDD5) close to this value with r2=0.357 (Supplementary Fig. 7) and an overall variation in CWM values explained by environment of 25.3% (cumulative variance of all 18 constrained axes in a RDA). The non-forest plots showed the same functional continua, but with lower total amount of variation in CWM accounted for by the first two PCA axes (41.8%, Supplementary Fig. 8) and much lower overall variation explained by environment. For non-forests, the best correlation of any CWM trait with environment was the one of volumetric content of coarse fragments in the soil (CoarseFrags) and leaf C content per dry mass with r2=0.042 (Supplementary Fig. 9). Similarly, the cumulative variance of all 18 constrained axes according to RDA was only 4.6%. This shows, on the one hand, that forest and non-forest vegetation are characterized by the same interrelationships of CWM traits, and on the other hand, that the relationships of CWM values with the environment were much stronger for forests than for non-forest formations. The coefficients of determination were even higher than those previously reported for trait-environment relationships for North American forests (between CWM of seed mass and maximum temperature, r2=0.281)3. 659 660 661 Resampling procedure in environmental space. In order to achieve a more even representation of plots across the global climate space, we first subjected the same 30 global climate and soil variables as described above, to a Principal Component Analysis (PCA), 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 using the climate space of the whole globe, irrespective of the presence of plots in this space, and scaling each variable to a mean of zero and a standard deviation of one. We used a 2.5 arc minute spatial grid, which comprised 8,384,404 terrestrial grid cells. We then counted the number of vegetation plots in the sPlot database that fell into each grid cell. For this analysis, we did not use the full set of 1,117,369 plots with trait information (see above), but only those plots that had a location inaccuracy of max. 3 km, resulting in a total of 799,400 plots. The resulting PCA scores based on the first two principal components (PC1-PC2) were rasterized to a 100 × 100 grid in PC1-PC2 environmental space, which was the most appropriate resolution according to a sensitivity analysis. This sensitivity analysis tested different grid resolutions, from a coarse-resolution bivariate space of 100 grid cells (10 × 10) to a very fineresolution space of 250,000 grid cells (500 × 500), iteratively increasing the number of cells along each principal component by 10 cells. For each iteration, we computed the total number of sPlot plots per environmental grid cell and plotted the median sampling effort (number of plots) across all grid cells versus the resolution of the PC1-PC2 space. We found that the curve flattens off at a bivariate environmental space of 100 × 100 grid cells, which was the resolution for which the median sampling effort stabilized at around 50 plots per grid cell. As a result, we resampled plots only in environmental cells with more than 50 plots (858 cells in total). 680 681 682 683 684 685 686 687 688 689 690 691 692 693 To optimize our resampling procedure within each grid cell, we used the heterogeneityconstrained random (HCR) resampling approach53. The HCR approach selects the subset of vegetation plots for which those plots are the most dissimilar in their species composition while avoiding selection of plots representing peculiar and rare communities that differ markedly from the main set of plant communities (outliers), thus providing a representative subset of plots from the resampled grid cell. We used the turnover component of the Jaccard’s dissimilarity index (βjtu54) as a measure of dissimilarity. The βjtu index accounts for species replacement without being influenced by differences in species richness. Thus, it reduces the effects of any imbalances that may exist between different plots due to species richness. We applied the HCR approach within a given grid cell by running 1,000 iterations of randomly selecting 50 plots out of the total number of plots available within that grid cell. Where the cell contained 50 or fewer plots, all were included and the resampling procedure was not run. This procedure thinned out over-sampled climate types, while retaining the full environmental gradient. 694 695 696 697 698 699 700 701 702 703 704 All 1,000 random draws of a given grid cell were subsequently sorted according to the decreasing mean of βjtu between pairs of vegetation plots and then sorted again according to the increasing variance in βjtu between pairs of vegetation plots. Ranks from both sortings were summed for each random draw, and the random draw with the lowest summed rank was considered as the most representative of the focal grid cell. Because of the randomized nature of the HCR approach, this resampling procedure was repeated 100 times for each of the 858 grid cells. This enabled us to produce 100 different subsamples out of the full sample of 799,400 vegetation plots subjected to the resampling procedure. Each of these 100 subsamples was finally subjected to ordinary linear regression, PCA and RDA as described above. We calculated the mean correlation coefficient across the 100 resampled data sets for each environmental variable with each trait. 705 706 707 708 709 710 711 712 713 714 715 To plot bivariate relationships, we used the mean intercept and slope of these relationships. PCA loadings of all 100 runs were stored and averaged. As different runs showed different orientation on the first PCA axes, we switched the signs of the axis loadings in some of the runs to make the 100 PCAs comparable to the reference PCA, based on the total data set. Across the 100 resampled data sets, we then calculated the minimum and maximum loading for each of the two PCA axes and plotted the result as ellipsoid. We also collected the posthoc regressions coefficients of PCA scores with the environmental variables in each of the 100 runs, switched the signs accordingly and plotted the correlations to PC1 and PC2 as ellipsoids. The result is a synthetic PCA of all 100 runs. To illustrate the coverage of plots in PCA space, we used plot scores of one of the 100 random runs. Similarly, the coefficients of determination obtained from the RDAs of these 100 resampled sets were averaged. 716 717 718 719 The mean PCA loadings across these 100 subsets (summarized in Supplementary Fig. 10) were fully consistent with those of the full data set in Fig. 2, with the same two functional continua in plant size and diaspore mass (from bottom left to top right), and perpendicular to that, the leaf economics spectrum. The variation in CWM accounted for by the first two axes was on average 50.9% ± 0.04 standard deviations (SD), and thus, virtually identical with that in the total dataset. In contrast, the variation explained on average by macroclimate and soil conditions (26.5% ± 0.01 SD as average cumulative variance of all 18 constrained axes in the RDAs across all 100 runs) was considerably larger than that for the total dataset, which is also reflected in consistently higher correlations between traits and environmental variables (Supplementary Fig. 11). The highest mean correlation was encountered for plant height and PET (mean r2=0.342 across 100 runs). PET was a better predictor for plant height than the precipitation of the wettest months (bio13, mean r2=0.231), as had been suggested previously6. The correlation of PET with stem specific density (mean r2=0.284) and warmth of the growing season (expressed as growing degree days above the threshold 5°C, GDD5) with leaf N:P ratio (mean r2=0.250) ranked among the best 12 correlations encountered out of all 540 trait-environment relationships, which confirms the patterns found in the whole data set (compared with Fig. 4). Overall, the coefficients of determination were much closer to the ones reported from other studies with a global collection of a few hundred plots (r2 values ranging from 36% to 53% based on multiple regressions of single traits with five to six environmental drivers22). 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 Data availability statement 738 739 740 741 742 The data contained in sPlot (the vegetation-plot data complemented by trait and environmental information) are available by request, through contacting any of the sPlot consortium members for submitting a paper proposal. The proposals should follow the Governance and Data Property Rules of the sPlot Working Group, which are available on the sPlot website (www.idiv.de/sPlot). 743 744 Acknowledgements 745 746 747 748 749 sPlot has been initiated by sDiv, the Synthesis Centre of the German Centre for Integrative Biodiversity Research (iDiv) Halle-Jena-Leipzig, funded by the German Research Foundation (FZT 118) and now is a platform of iDiv. H.B., J.De., O.Pu, U.J., B.J.-A., J.K., D.C., F.M.S., M.W. and C.W. appreciate direct funding through iDiv. For all further acknowledgements see the Supplementary Information. 750 751 Author contributions 752 753 754 755 756 757 758 759 760 761 762 763 H.B. and U.J. wrote the first draft of the manuscript, with considerable input by B.J.-A. and R.F.; H.B. carried out most of the statistical analyses and produced the graphs; H.B., O.Pu. and U.J. initiated sPlot as an sDiv working group and iDiv platform; J.De. compiled the plot databases globally; J.De., S.M.H., U.J., O.Pu. and F.J. harmonized vegetation databases; J.De. and B.J.-A. coordinated the sPlot consortium; J.K. provided the trait data from TRY; F.S. performed the trait data gap filling; O.Pu. produced the taxonomic backbone; B.J.-A., G.S. and E. Welk compiled environmental data and produced the global maps; S.M.H. wrote the Turboveg v3 software, which holds the sPlot database; J.L. and T.H. wrote the resampling algorithm. Many authors participated in one or more of the three sPlot workshops at iDiv where the sPlot initiative was conceived and planned, and evaluation of the data and first drafts were discussed. All other authors contributed data. All authors contributed to writing the manuscript. 764 765 Declaration of competing interests 766 The authors declare no competing interests. 767 768 References 769 770 1. Warming, E. Lehrbuch der ökologischen Pflanzengeographie – Eine Einführung in die Kenntnis der Pflanzenvereine. (Borntraeger, Berlin, 1896). 771 772 2. Garnier, E. et al. Plant functional markers capture ecosystem properties during secondary succession. Ecology 85, 2630–2637 (2004). 773 774 3. Ordoñez, J.C. et al. A global study of relationships between leaf traits, climate and soil measures of nutrient fertility. Global Ecol. Biogeogr. 18, 137–149 (2009). 775 776 4. Garnier, E., Navas, M.-L. & Grigulis, K. Plant functional diversity - Organism traits, community structure, and ecosystem properties. (Oxford Univ. Press, 2016). 777 5. Díaz, S. et al. The global spectrum of plant form and function. Nature 529, 167-171 (2016). 778 6. Moles, A.T. et al. Global patterns in plant height. J. Ecol. 97, 923-932 (2009). 779 7. Wright, I.J. et al. The worldwide leaf economics spectrum. Nature 428, 821-827 (2004). 780 781 8. Reich, P.B. The world–wide ‘fast–slow’ plant economics spectrum: a traits manifesto. J. Ecol. 102, 275–301 (2014). 782 783 9. Adler, P.B. et al. Functional traits explain variation in plant life history strategies. Proc. Natl. Acad. Sci. USA 111, 740–745 (2014). 784 785 10. Marks, C.O. & Lechowicz, M.J. Alternative designs and the evolution of functional diversity. Am. Nat. 167, 55–67 (2006). 786 787 11. Grime, J.P. Trait convergence and trait divergence in herbaceous plant communities: Mechanisms and consequences. J. Veg. Sci. 17, 255–260 (2006). 788 789 12. Muscarella, R. & Uriarte, M. Do community-weighted mean functional traits reflect optimal strategies? Proc. R. Soc. B. 283, 20152434 (2016). 790 791 13. Swenson, N.G. & Weiser, M.D. Plant geography upon the basis of functional traits: An example from eastern North American trees. Ecology 91, 2234–2241 (2010). 792 793 14. Fyllas, N.M. et al. Basin-wide variations in foliar properties of Amazonian forest: phylogeny, soils and climate. Biogeosciences 6, 2677-2708 (2009). 794 795 15. Swenson, N.G. et al. Phylogeny and the prediction of tree functional diversity across novel continental settings. Global Ecol. Biogeogr. 26, 553–562 (2017). 796 797 16. Swenson, N.G. et al. The biogeography and filtering of woody plant functional diversity in North and South America. Global Ecol. Biogeogr. 21, 798–808 (2012). 798 17. Wright, I.J. et al. Global climatic drivers of leaf size. Science 357: 917–921 (2017). 799 800 18. Mayfield, M.M. & Levine, J.M. Opposing effects of competitive exclusion on the phylogenetic structure of communities. Ecol. Lett. 13, 1085 – 1093 (2010). 801 802 19. Kraft, N.J.B. et al. Community assembly, coexistence and the environmental filtering metaphor. Funct. Ecol. 29, 592–599 (2015). 803 804 20. Barboni, D. et al. Relationships between plant traits and climate in the Mediterranean region: A pollen data analysis. J. Veg. Sci. 15, 635-646 (2004). 805 806 21. Borgy, B. et al. Plant community structure and nitrogen inputs modulate the climate signal on leaf traits. Global Ecol. Biogeogr. 26, 1138-1152 (2017). 807 808 809 22. van Bodegom, P.M, Douma, J.C. & Verheijen, L.M. A fully traits-based approach to modeling global vegetation distribution. Proc. Natl. Acad. Sci. USA 111, 13733–13738 (2014). 810 811 23. Moles, A.T. et al. Which is a better predictor of plant traits: Temperature or precipitation? J. Veg. Sci. 25, 1167–1180 (2014). 812 813 24. Ordoñez, J.C. et al. Plant strategies in relation to resource supply in mesic to wet environments: Does theory mirror nature? Am. Nat. 175, 225–239 (2010). 814 815 816 25. Simpson, A.J., Richardson, S.J. & Laughlin, D.C. Soil–climate interactions explain variation in foliar, stem, root and reproductive traits across temperate forests. Global Ecol. Biogeogr. 25, 964-978 (2016). 817 818 26. Lienin, P. & Kleyer, M. Plant leaf economics and reproductive investment are responsive to gradients of land use intensity. Agric. Ecosyst. Environ. 145, 67-76 (2011). 819 820 821 27. Maire, V. et al. Habitat filtering and niche differentiation jointly explain species relative abundance within grassland communities along fertility and disturbance gradients. New Phytol. 196, 497–509 (2012). 822 823 824 28. Craine, J.M. et al. Global patterns of foliar nitrogen isotopes and their relationships with climate, mycorrhizal fungi, foliar nutrient concentrations, and nitrogen availability New Phytol. 183, 980–992 (2009). 825 826 29. Güsewell, S. N:P ratios in terrestrial plants: variation and functional significance. New Phytol. 164, 243–266 (2004). 827 828 30. Reich, P.B. & Oleksyn, J. Global patterns of plant leaf N and P in relation to temperature and latitude, Proc. Natl. Acad. Sci. USA 101, 11001-11006 (2004). 829 830 31. Scheiter, S., Langan, L. & Higgins, S.I. Next generation dynamic global vegetation models: learning from community ecology. New Phytol. 198, 957-969 (2013). 831 832 32. Boyle, B. et al. The Taxonomic Name Resolution Service: an online tool for automated standardization of plant names. BMC Bioinformatics 14, 16 (2013). 833 834 33. Bremer, B. et al. An update of the Angiosperm Phylogeny Group classification for the orders and families of flowering plants: APG III. Bot. J. Linn. Soc. 161, 105–121 (2009). 835 836 837 34. Schrodt, F. et al. BHPMF – a hierarchical Bayesian approach to gap-filling and trait prediction for macroecology and functional biogeography. Global Ecol. Biogeogr. 24, 1510– 1521 (2015). 838 839 35. Kattge J. et al. TRY—a global database of plant traits. Glob. Change Biol. 17, 2905–2935 (2011). 840 841 842 36. Shan, H. et al. Gap filling in the plant kingdom - Trait prediction using hierarchical probabilistic matrix factorization. Proceedings of the 29th International Conference for Machine Learning (ICML 2012) 1303−1310 (2012). 843 844 845 37. Fazayeli, F. et al. Uncertainty quantified matrix completion using Bayesian Hierarchical Matrix factorization. 13th International Conference on Machine Learning and Applications (ICMLA 2014) 312−317 (2014). 846 847 848 38. Borgy, B. et al. Sensitivity of community-level trait–environment relationships to data representativeness: A test for functional biogeography. Global Ecol. Biogeogr. 26, 729–739 (2017). 849 850 39. Herz, K. et al. Drivers of intraspecific trait variation of grass and forb species in German meadows and pastures. J. Veg. Sci. 28, 705–716 (2017). 851 852 40. Karger, D.N. et al. Climatologies at high resolution for the earth’s land surface areas. Sci. Data 4, 170122. doi: 10.1038/sdata.2017.122 (2017) 853 854 855 41. Karger, D.N. et al. Climatologies at high resolution for the Earth land surface areas (Version 1.1). World Data Center for Climate (WDCC) at DKRZ. http://chelsaclimate.org/downloads/ (2016). 856 857 42. Dee, D. P. et al. The ERA-Interim reanalysis: configuration and performance of the data assimilation system. Q.J.R. Meteorol. Soc. 137, 553–597 (2011). 858 859 860 43. Synes, N.W. & Osborne, P.E. Choice of predictor variables as a source of uncertainty in continental-scale species distribution modelling under climate change. Global Ecol. Biogeogr. 20, 904–914 (2011). 861 862 863 44. Enquist, B. et al. Scaling from traits to ecosystems: developing a general trait driver theory via integrating trait-based and metabolic scaling theories. Adv. Ecol. Res. 52, 249-318 (2015). 864 865 45. Buzzard, V. et al. Re-growing a tropical dry forest: functional plant trait composition and community assembly during succession. Funct. Ecol. 30, 1006–1013 (2016). 866 867 46. Rao, C.R. Diversity and dissimilarity coefficients: A unified approach. Theor. Popul. Biol. 21, 24–43 (1982). 868 869 47. Champely, S. & Chessel, D. Measuring biological diversity using Euclidean metrics. Environ. Ecol. Stat. 9, 167-177 (2002). 870 871 48. Dowle, M. et al. data.table: Extension of Data.frame. R package version 1.9.6. (2015) https://CRAN.R-project.org/package=data.table 872 873 874 49. Hawkins, B.A. et al. Structural bias in aggregated species-level variables driven by repeated species co-occurrences: a pervasive problem in community and assemblage data. J. Biogeogr. 44, 1199–1211 (2017). 875 876 50. Knijnenburg, T.A. et al. Fewer permutations, more accurate P-values. Bioinformatics 25, i161–i168 (2009). 877 878 51. Friendly, M. Corrgrams: Exploratory displays for correlation matrices. Am. Statistician, 56, 316–324 (2002). 879 880 52. Oksanen, J. et al. vegan: Community Ecology Package. R package version 2.3-3 (2016). https://CRAN.R-project.org/package=vegan 881 882 53. Lengyel, A., Chytrý, M. & Tichý, L. Heterogeneity-constrained random resampling of phytosociological databases. J. Veg. Sci. 22, 175–183 (2011). 883 884 54. Baselga, A. The relationship between species replacement, dissimilarity derived from nestedness, and nestedness. Global Ecol. Biogeogr. 21, 1223-1232 (2012). 885 886 55. Garnier, E. et al. Towards a thesaurus of plant characteristics: an ecological contribution. J. Ecol. 105, 298-309 (2017). www.top-thesaurus.org 887 888 889 Tables 890 891 892 893 894 895 896 Table 1: Traits used in this study and their function in the community. Traits are arranged according to the degree to which they should respond to macroclimatic drivers. ↑↓ in the trait column denotes opposing relationships,  in the description column denotes trade-offs. For trait units, plot-level trait means and within-plot trait variance see Table 2. Trait Description Function Expected correlation with macroclimate Specific leaf area, Leaf area, Leaf fresh mass, Leaf N, Leaf P ↑↓ Leaf dry matter content, Leaf N per area, Leaf C Leaf economics spectrum7-8,17: Thin, N-rich leaves with high turnover and high mass-based assimilation rates  Thick, N-conservative, long-lived leaves with low mass-based assimilation rates Productivity, competitive ability Very high12- Stem specific density Fast growth  Mechanical support, Longevity Productivity, drought tolerance Very high12,22 Conduit element length ↑↓ Stem conduit density Efficient water transport  Safe water transport Water use efficiency High Plant height Mean individual height of adult plants Competitive ability High6,12 Seed number per reproductive unit ↑↓ Seed mass, Seed length, Dispersal unit length Seed economics spectrum23: Small, well dispersed seeds  Seeds with storage reserve to facilitate establishment and increase survival Dispersal, regeneration Moderate23-24 Leaf N:P ratio P limitation (N:P > 15) N limitation (N:P < 10)29 Nutrient supply Moderate30 Leaf nitrogen isotope ratio (leaf δ15N) Access to N derived from N2 fixation  N supply via mycorrhiza Nitrogen source, soil depth Moderate28 13,17,21,23 897 898 899 900 901 902 903 904 905 906 907 Table 2: Traits, abbreviation of trait names, identifier in the Thesaurus Of Plant characteristics (TOP)55, units of measurement, observed values (obs.) standardized effect sizes (SES) and significance (p) of SES for means and variances of both plot-level trait means (community-weighted means, CWMs) and within-plot trait variances (community-weighted variances, CWVs). CWMs and CWVs were based on gap-filled traits for 1,115,785 and 1,099,463 plots, respectively. All trait values were loge-transformed prior to analysis and observed values are on the loge scale. SES are also based on loge-transformed values. Stem specific density is stem dry mass per stem fresh volume, specific leaf area is leaf area per leaf dry mass, leaf C, N and P are leaf carbon, nitrogen and phosphorus content, respectively, per leaf dry mass, leaf dry matter content is leaf dry mass per leaf fresh mass, leaf delta 15N is the leaf nitrogen isotope ratio, stem conduit density is the number of vessels and tracheids per unit area in a cross section, conduit element length refers to both vessels and tracheids. SESs were calculated by randomizing trait values across all species globally 100 times and calculating CWM and CWV with random trait values, but keeping all species abundances in plots (see Fig. 1). Tests for significance of SES were obtained by fitting generalized Pareto-distribution of the most extreme random values and then estimating p values form this fitted distribution50. * indicates significance at p < 0.05. CWM Trait Leaf area Specific leaf area Leaf fresh mass Leaf dry matter content Leaf C Leaf N Leaf P Leaf N per area Leaf N:P ratio Leaf δ15N Seed mass Seed length Seed number per reproductive unit Dispersal unit length Abbreviation LA SLA Leaf.fresh.mass LDMC LeafC LeafN LeafP LeafN.per.area Leaf.N:P.ratio Leaf.delta15N Seed.mass Seed.length Seed.num.rep.unit TOP 25 50 35 45 452 462 463 481 103 91 - Unit mm2 m2 kg-1 g g g-1 mg g-1 mg g-1 mg g-1 g m-2 g g-1 ppm mg mm obs. 6.130 2.850 -2.125 -1.294 6.116 3.038 0.535 0.251 2.444 0.521 0.407 1.069 6.179 Disp.unit.length 90 mm 1.225 mean SES CWV -9.75 9.89 -13.28 -5.67 -3.77 4.22 9.57 -9.06 -11.95 -3.58 -11.19 -4.51 7.67 p * * * * * * * * * * * * * variance obs. SES p 1.691 12.53 * 0.172 12.88 * 1.395 10.83 * 0.101 11.52 * 0.003 8.80 * 0.055 6.29 * 0.097 2.81 * 0.075 8.18 * 0.040 0.40 n.s. 0.254 6.68 * 2.987 3.69 * 0.294 5.50 * 2.783 4.40 * -2.51 * 0.343 6.50 * mean SES obs. 1.565 0.150 1.520 0.130 0.002 0.063 0.117 0.099 0.081 0.455 2.784 0.365 5.156 0.451 -3.21 -2.59 -1.33 -2.05 0.95 -1.78 -3.19 -5.17 -0.28 -2.74 2.82 -9.06 -4.67 1.44 p * n.s. * n.s. * * * n.s. * * * * n.s. variance obs. SES p 2.448 -0.27 n.s. 0.023 1.10 n.s. 2.311 0.01 n.s. 0.017 6.73 * 0.000 -0.38 n.s. 0.004 -0.13 n.s. 0.014 -2.11 * 0.010 1.54 n.s. 0.007 -0.39 n.s. 0.207 2.44 * 7.750 -2.81 * 0.134 -3.07 * 26.588 2.25 * * 0.203 -1.39 n.s. Plant height Plant.height Stem specific density SSD Stem conduit density Stem.cond.dens Conduit element length Cond.elem.length Mean SES Mean absolute SES 908 909 68 286 - m g cm-3 mm-2 µm -0.315 -0.869 4.407 5.946 -12.15 -14.93 15.08 -7.09 -3.50 8.66 * * * * 1.532 0.041 0.656 0.182 13.34 13.15 8.45 9.14 8.06 8.06 * * * * 1.259 0.058 0.975 0.367 -9.01 2.09 -0.95 7.12 -1.76 3.36 * * n.s. * 1.585 0.003 0.951 0.135 9.68 2.99 1.10 5.29 1.25 2.43 * * n.s. * 910 Captions of Figures 911 912 913 914 915 916 917 918 919 920 921 922 923 Fig. 1: Conceptual figure to illustrate Hypothesis 1, stating that environmental or biotic filtering of community trait values result in a) higher than expected variation of communityweighted means and b) lower than expected community-weighted variances of trait values. Both figures give an example for a single trait and show the relative abundance of trait values of all species in a plot. Black curves refer to observed plot-level trait values in two exemplary plots, while grey curves show plot-level trait values obtained from randomizing trait values across all species globally (see Methods). Randomization was done 100 times, but only one randomization event is shown. Deviation from random expectation was assessed with standardized effect sizes (SESs) for a) the variance in CWMs and b) for the mean in CWVs. Evidence for filtering is given in a) if the variance in plot-level trait means was higher than expected by chance (SES significantly positive) or b) if within-plot trait variance was typically lower than expected by chance (SES significantly negative, see Methods). 924 925 926 927 928 929 930 931 932 933 934 935 936 Fig. 2: Principal Component Analysis of global plot-level trait means (community-weighted means, CWMs). The plots (n=1,114,304) are shown by coloured dots, with shading indicating plot density on a logarithmic scale, ranging from yellow with 1–4 plots at the same position to dark red with 251–1142 plots. Prominent spikes are caused by a strong representation of communities with extreme trait values, such as heathlands with ericoid species with small leaf area and seed mass. Post-hoc correlations of PCA axes with climate and soil variables are shown in blue and magenta, respectively. Arrows are enlarged in scale to fit the size of the graph; thus, their lengths show only differences in variance explained relative to each other. Variance in CWM explained by the first and second axis was 29.7% and 20.1%, respectively. The vegetation sketches schematically illustrate the size continuum (short vs. tall) and the leaf economics continuum (low vs. high LDMC and leaf N content per area in light and dark green colours, respectively). See Table 2 and Supplementary Table 2 for the description of traits and environmental variables. 937 938 939 940 941 942 943 944 945 946 947 948 949 950 Fig. 3: Principal Component Analysis of global within-plot trait variances (communityweighted variances, CWVs). The plots (n=1,098,015) are shown by coloured dots, with shading indicating plot density on a logarithmic scale, ranging from yellow with 1–2 plots at the same position to dark red with 631–1281 plots. Post-hoc correlations of PCA axes with climate and soil variables are shown in blue and magenta, respectively. Arrows are enlarged in scale to fit the size of the graph; thus, their lengths show only differences in variance explained relative to each other. Variance in CWV explained by the first and second axis was 24.9% and 13.4%, respectively. CWV values of all traits increased from the left to the right, which reflects increasing species richness (r2 = 0.116 between scores of the first axis and number of species in the communities for which traits were available). The vegetation sketches schematically illustrate low and high variation in the plant size and leaf economics continua. See Table 2 and Supplementary Table 2 for the description of traits and environmental variables. 951 952 953 954 955 956 957 958 959 960 961 962 963 Fig. 4: The two strongest relationships found for global plot-level trait means (communityweighted means, CWMs) in the sPlot dataset. CWM of the natural logarithm of stem specific density [g cm-3] as a) global map, interpolated by kriging within a radius of 50 km around the plots using a grid cell of 10 km, and b) function of potential evapotranspiration (PET, r2=0.156). CWM of the natural logarithm of the N:P ratio [g g-1] as c) global kriging map and d) function of the warmth of the growing season, expressed as growing degree days over a threshold of 5°C (GDD5, r2=0.115). Plots with N:P ratios > 15 (of 2.71 on the loge scale) tend to indicate phosphorus limitation29 and are shown above the broken line in red colour (90,979 plots, 8.16% of all plots). The proportion of plots with N:P ratios > 15 increases with GDD5 (r2=0.895 for a linear model on the log response ratio of counts of plots with N:P > 15 and ≤15 counted within bins of 500 GDD5). a) random Relative abundance observed ( . = ) − . . Trait value b) Relative abundance ( observed random Trait value . ) = − . . a b CWM Stem specific density [g cm -3 ] 0.8 0.2 c d CWM Leaf N:P ratio [g g -1 ] 29 6