Plant functional and taxonomic diversity in European grasslands along climatic gradients

1Department of Environmental Science, Institute for Water and Wetland Research, Radboud University, Nijmegen, The Netherlands 2Department of Biology and Biotechnologies “Charles Darwin”, Sapienza Università di Roma, Rome, Italy 3National Research Council, Institute of Research on Terrestrial Ecosystems (CNRIRET), Monterotondo, Italy 4Department of Aquatic Ecology and Environmental Biology, Institute for Water and Wetland Research, Radboud University, Nijmegen, The Netherlands 5Vegetation Ecology Group, Institute of Natural Resource Sciences (IUNR), Zurich University of Applied Sciences (ZHAW), Wädenswil, Switzerland 6Plant Ecology, Bayreuth Center of Ecology and Environmental Research (BayCEER), Bayreuth, Germany 7German Centre for Integrative Biodiversity Research (iDiv) HalleJenaLeipzig, Leipzig, Germany 8WSL Swiss Federal Research Institute, Birmensdorf, Switzerland 9Plant Biology and Ecology, Faculty of Science and Technology, University of the Basque Country UPV/EHU, Bilbao, Spain 10EnvixLab, University of Molise, Pesche (IS), Italy 11Department of Biotechnologies and Life Sciences (DBSV), University of Insubria, Varese, Italy 12Department of Botany and Zoology, Faculty of Science, Masaryk University, Brno, Czech Republic 13Institute of Biology, Martin Luther University HalleWittenberg, Halle (Saale), Germany 14German Centre for Integrative Biodiversity Research (iDiv) HalleJena Leipzig, Leipzig, Germany 15General Vegetation Laboratory, Komarov Botanical Institute RAS, SaintPetersburg, Russia 16Phytodiversity Problems Laboratory, Samara Federal Research Center RAS, Institute of the Ecology of the Volga Basin RAS, Togliatti, Russia 17EnvixLab, University of Molise, Termoli, Italy 18School of Environment, Earth and Ecosystem Sciences, Faculty of Science, Technology, Engineering and Mathematics, The Open University, Milton Keynes, UK 19Faculty of Geography and Earth Sciences, University of Latvia, Rīga, Latvia


| INTRODUC TI ON
Grasslands are among the most diverse ecosystems in Europe (Wilson et al., 2012;Dengler et al., 2020). They have been extensively studied for a long time and with long-term monitoring schemes (Scholz, 1975;Willems, 1983;Tilman et al., 2006). Various assembly processes have been put forward that may explain the origin and maintenance of European grassland diversity, e.g., competitive hierarchy and niche partitioning (Mori et al., 2018). Yet, over the last 50 years, European grassland diversity has seen a dramatic decrease, often being attributed to increased nutrient availability (Wesche et al., 2012), overgrazing (Dengler et al., 2020) or drought (Carmona et al., 2012;Nogueira et al., 2018). Insight into patterns of plant diversity over environmental gradients is needed to aid deeper understanding of the effects of global change on biodiversity (e.g., Mooney et al., 2009;Cardinale et al., 2012;Funk et al., 2017) and may also improve our understanding of the mechanisms that underlie community assembly (MacArthur & Levins, 1967;McGill et al., 2006).
In the context of macroecology, diversity patterns of vegetation are typically discussed via various filtering mechanisms (e.g., Weiher et al., 2011). First, the dispersal filter determines the ability of a species to be present in a specific location, and hence the regional plant species pool and plant trait pool (Cadotte & Tucker, 2017). Second, environmental conditions act as an additional filter on plant communities, sorting those that fulfill local (fundamental) niche requirements constituted by physiological constraints (Leibold et al., 2004;Tingley et al., 2014). Under the favorability hypothesis, the more extreme or unfavorable environmental conditions are, the more selective environmental filters are (Fischer, 1960). This suggests that only plants with trait values well adapted to the extreme Funding information CCFB, MAJH and SH were supported by the ERC project (CoG SIZE 647224), IB by the Basque Government (IT936-16), MC by the Czech Science Foundation (19-28491X), and SR by the University of Latvia grant (AAp2016/B041//Zd2016/AZ03).

Co-ordinating Editor: Meelis Pärtel
Methods: We quantified functional and taxonomic richness of 55,748 vegetation plots. Six plant traits, related to resource acquisition and conservation, were analysed to describe plant community functional composition. Using a null-model approach we derived functional richness effect sizes that indicate higher or lower diversity than expected given the taxonomic richness. We assessed the variation in absolute functional and taxonomic richness and in functional richness effect sizes along gradients of minimum temperature, temperature range, annual precipitation, and precipitation seasonality using a multiple general additive modelling approach. conditions persist, resulting in a reduction in community functional diversity in extreme conditions (de Bello et al., 2009;Mayfield & Levine, 2010;Shen et al., 2016). A third filter may be biotic interactions, representing a counter-gradient with competitive interspecific interactions being the main driver of plant community composition at low abiotic stress (the stress gradient hypothesis; Bertness & Callaway, 1994). While this competition may lead to exclusion of species and reduced trait variation (Grime, 2006;Mayfield & Levine, 2010;Kunstler et al., 2012), facilitative interactions have been suggested to dominate when abiotic stress is high (Bertness & Callaway, 1994;Brooker & Callaghan, 1998), leading to higher trait divergence and increased functional diversity (Valiente-Banuet & Verdú, 2007;McIntire & Fajardo, 2014). Finally, the occurrence and abundance of a plant may also be influenced by temporal variation in climate or other stressors (Díaz & Cabido, 1997;González-Moreno et al., 2015;Fischer et al., 2020). Plants can occur in a place as long as species' niche requirements are met at some time of the year, i.e., exploiting temporarily empty niches (Godoy et al., 2009). This expansion of the realized niche can be especially important for the diversity of plant communities in temperate regions (Scheiner & Rey-Benayas, 1994;Breitschwerdt et al., 2018).

Results
When describing diversity patterns, the focus should not only be on plant species identity (i.e., their taxonomy) but also on plant traits. As traits describe a more direct link between the performance of an organism, local environmental conditions, and a plant's functioning in the community (Keddy, 1992;Funk et al., 2017), they play a critical role in determining local plant community diversity (Tilman et al., 1997;Weiher et al., 1998). Especially when non-random processes determine community assemblages, species diversity is not an adequate surrogate for functional diversity (Díaz & Cabido, 2001). However, both species and trait composition vary along environmental gradients (e.g., Raymundo et al., 2018), making the separation of taxonomic and functional diversity difficult (e.g., Hooper et al., 2002;Petchey & Gaston, 2002). To control for potentially coincidental similarity in diversity-environment trends, specifically designed field experiments are needed or one can use null models that check for functional diversity patterns different from random sampling of species (Gotelli & Graves, 1996;Swenson et al., 2012).
Taken together, the net influence of the various filtering mechanisms on taxonomic and functional plant richness along large environmental stress gradients remains elusive. Part of this ambiguous question is the value of considering plant trait diversity in addition to taxonomic diversity. Here, we analyzed functional and taxonomic richness in European grasslands along temperature and precipitation gradients. We quantified the absolute functional and taxonomic richness directly derived from the community data, as well as the effect size of functional richness. The effect size indicates higher or lower functional richness than expected given the observed taxonomic richness (e.g., Harvey et al., 1983;Götzenberger et al., 2016).
First, we focused on annual minimum temperature, as extreme temperatures limit metabolic activity, growth and constrain leaf size, and thus are likely to filter trait composition (Went, 1953). Second, we included annual precipitation, as drought constrains leaf longevity, photosynthetic efficiency, and seed mass (Sandel et al., 2010). Last, we included annual temperature range and precipitation seasonality to be able to assess the influence of seasonality on plant diversity at the continental scale (Scheiner & Rey-Benayas, 1994;Godoy et al., 2009  . All grassland habitats (EUNIS group R) and coastal dune habitats (EUNIS habitats N15, N16 and N17) were selected. Plots not classified to any of these habitat types but assigned in the EVA database to the phytosociological class Molinio-Arrhenatheretea (Mucina et al., 2016) were added to the data set.
Data covered the whole of Europe. The islands of Macaronesia were excluded to remove any effects of isolated oceanic islands with significant within-island speciation and endemism (Humphries, 1979). We selected plots recorded between 1979 and 2013 to create an optimal match with climatic data. In addition, we selected plots with a georeferencing uncertainty smaller than 1 km, which further improved the match with high-resolution climatic data. Such a direct link between plant community and climate data reduces the effects of confounding factors like habitat heterogeneity (Szilágyi & Meszéna, 2009). We only included flowering plants as the available traits for this group are consistent within the data set. By doing so, we excluded ferns, bryophytes, lichens and fungi.
These selection criteria resulted in a georeferenced occurrence data set with presences of all selected plant species in all selected plots. From this initial data set, we excluded plots with particular ecological conditions (Münkemüller et al., 2020) using the following criteria: (a) plots influenced by high salinity, identified by the presence of species known to favor extreme salty conditions (e.g., Suaeda maritima) were excluded, as their composition is likely driven by salinity rather than temperature or precipitation; and (b) plots with a cover of tree or tall shrub species, identified as all plant species with an average height of more than 3 m (obtained from the TRY database, see below), that exceeded 5% in cover. All sources of vegetation-plot data are listed in Appendix S1.

| Trait data
We obtained trait data from the TRY database version 5 (Kattge et al., 2020), containing both public and restricted data sets. We retrieved root trait data from the Fine-Root Ecology Database (FRED; Iversen et al., 2018). All data sources are listed in Appendix S2. After standardization of trait units, we merged trait data from the two databases and calculated species averages. While averaging excludes intraspecific trait variation and trait plasticity, which determines a large part of community diversity (Albert et al., 2010;Jung et al., 2010;Ross et al., 2017;Barbour et al., 2019;Niu et al., 2020), this step was necessary as trait data were not available for most species across the different plots.
Functional richness depends on the selected set of traits. We used a sequence of selection criteria leading to our final trait data set. First, we selected traits based on their importance for vegetative growth and reproduction (Wright et al., 2004;Díaz et al., 2016).
Then, we selected both above-and below-ground traits with good data coverage. Third, we removed traits of the same plant organ with high covariance (e.g., specific leaf area and leaf dry matter content) in order to equalize organ and trait importance in the functional diversity calculation, keeping the trait with the highest data coverage.
Fourth, we selected numeric traits only to minimize the effects of dimension choice to calculate functional diversity (Maire et al., 2015).
Our final trait data set contains six plant traits (Table 1).
As functional richness calculations (see below) are sensitive to the completeness of the species list per plot (Pakeman, 2014) and thus require a complete trait data set, we filled the gaps (14.7% of overall missing data) in our trait database (Table 1; Appendix S3).
Imputing missing data is considered a better alternative to removing incomplete plots, which can potentially result in systematic biases (Nakagawa & Freckleton, 2008). We used the mice function of the mice package in R-3.6.1, which relies on between-trait correlations using multi-variate imputation with chained equations to fill gaps. This method has been shown to have lower imputation error and bias compared to other multiple imputation approaches, especially when the percentage of missing data is high (Penone et al., 2014). Although gap-filling techniques require only one trait value per species, we only retained species with at least two trait values to achieve a more accurate imputation. This resulted in the removal of complete relevés if the trait values for at least one species in the community did not meet this requirement. This resulted in the data set having 61,714 relevés containing a total of 2,884 species. Prior to the imputation, we log-transformed all trait data since some approaches could be sensitive to data with varying scales in the variables (Penone et al., 2014). Then, we used the predictive mean matching method for the imputation of all traits to preserve nonlinear trait-trait relationships. We ran 24 multiple imputations (van Buuren, 2012) where trait predictions were updated five times in the chained equations to improve the quality of the predictions.
Because multiple imputation performance benefits from the addition of other traits that may be related to those of interest, we also included plant longevity (categorical trait), as this trait was complete at the genus level. The gap-filling step resulted in 24 trait data sets, where the variation between them represents the uncertainty of the gap-filling technique. Qualitatively consistent results from different imputed data sets indicate that the conclusions on the directionality of the effects are not influenced by the gap-filling procedure.

| Climate data
We constructed four environmental gradients. The first was daily average minimum temperature (°C) of the coldest month of the year and the second was annual precipitation (mm). For both, environments were assumed to be harsher at the low end of the gradient (Went, 1953). Third, we considered precipitation seasonality (%), defined as the coefficient of variation based on monthly precipitation. Fourth, we considered annual temperature range (°C), defined as the absolute difference between the daily average maximum temperature of the warmest month and the daily average minimum TA B L E 1 List of functional plant traits used. Missing data for species represent the percentage of species that have no data for that trait. Missing data overall represent the percentage of gaps in the overall data set for that trait

| Diversity
We used two diversity metrics, taxonomic richness and multi-variate number of species per plot occurred along the entire environmental gradient and followed the same pattern as plots with more than six species (Appendix S4). Removing these communities resulted in 24 final data sets, each with 55,748 plots and 2,830 different species ( Figure 1).

| Null model
We used a null model to test for potential effects of taxonomic richness on functional richness, revealing non-random patterns that indicate ecological processes as opposed to random patterns expected by chance (Gotelli & Graves, 1996;Swenson et al., 2012;Chalmandrier et al., 2013). Hence, we calculated a functional richness effect size using results from a null model which randomly shuffles species trait values in the trait data set but retains species' trait combinations. Specifically, the species' PCoA scores were randomly shuffled among species (but not within species) from the entire species pool .
The number of species per community was kept equal to the number of species observed to ensure the assessment of functional richness trends, independent of differences in taxonomic richness (similarly to Swenson et al., 2012;Craven et al., 2018).
The effect sizes were calculated using probabilities since the dis-

| Statistical analyses
We used the average of the original functional richness values and effect size values in one grid cell (~1 km 2 ) if more than one plot was sampled to reduce possible effects of spatial autocorrelation. This resulted in 19,179 data points and grid cells (Figure 1). For each of the 24 imputed data sets, functional richness, taxonomic richness, and effect sizes were assessed along gradients of minimum temperature, temperature range, annual precipitation and precipitation seasonality in multiple generalized additive models with a Gaussian distribution using the gamm4 package. These models allowed for finding patterns in the data without a priori hypotheses on the shape of the relationship. However, we restricted the curviness of the trends by setting the knots (k) to 4 to prevent overfitting. Since the difference between the models was small (Appendix S6), we averaged the trait data over the 24 imputed data sets, recalculated functional richness values, reran the null model, and used these data in our plots of functional richness and effect sizes so that confidence intervals could indicate model uncertainty instead of differences between imputed data sets.

| RE SULTS
The diversity of grassland communities varied considerably over temperature and precipitation gradients (Figure 2; Appendix S7). We found functional and taxonomic richness to be lower in areas with low minimum temperatures and in areas with narrow temperature ranges (Figure 2a,c). These corresponded to locations in mountain ranges or with higher latitude or longitude and locations on islands or along the west coast of France, Belgium and the Netherlands, respectively (Appendix S9). This similarity in trend between functional and taxonomic richness may be due to a substantial correlation (r = 0.7; Appendix S8). However, at wide annual temperature ranges, functional richness increased with increasing minimum temperatures (Figure 2a), while taxonomic richness decreased with increasing minimum temperatures after reaching an optimum at a minimum temperature of around −5 °C (Figure 2c). These plots were located in continental Europe, where the highest minimum temperatures (i.e., highest functional richness values) were found along the coasts of southern Europe (Appendix S9). Note that for all temperature variables (annual maximum, maximum of the warmest quarter, annual minimum, minimum of the coldest quarter, and annual mean temperature), both richness metrics showed the same non-linear pattern as for the minimum temperature (Appendix S10).
Likewise, temperature variation variables (annual temperature range and temperature seasonality) showed the same non-linear pattern (Appendix S10). Considering the variation of functional and taxonomic richness across precipitation gradients, we found a decrease In contrast, effect sizes did not show a clear pattern along the temperature range gradient (Figure 2e), nor did they vary much across both precipitation gradients (Figure 2f).

| D ISCUSS I ON
The results of our continental-scale study show that European grassland communities vary substantially in functional and taxonomic richness along temperature and precipitation gradients.
Highest richness in both traits and species was found in communi-  (Fischer, 1960).
Additionally, this pattern is in line with taxonomic and functional diversity patterns found in non-European grasslands (Moradi & Oldeland, 2019). Note that Moradi & Oldeland (2019) reported that precipitation limits plant diversity more than minimum temperature in dry areas, which contrasts with our results (Appendix S11).
This may suggest that, when drier areas would have been included in this study, taxonomic richness may decrease more toward more water-limited sites and determine the diversity patterns more than temperature. Another explanation of the peak in functional and taxonomic richness at average minimum temperature values is the mid-domain effect (Colwell & Hurtt, 1994;Colwell & Lees, 2000).
Here, the geometric limits in relation to the distribution of species increases the overlap of species ranges and thus of species richness in the middle of the gradient (Colwell & Lees, 2000). We assume this effect to be of minor importance as grasslands extend from southwestern Europe to Middle and Central Asia, many grassland species have large geographical ranges (extending beyond the area included in this analysis), and many included species also occur in other habitats besides grasslands.
Remarkably, functional richness values further increased toward higher minimum temperatures and wider temperature ranges, while taxonomic richness decreased. Environments with a wide temperature range may hold co-existing species with dissimilar trait values as they can exploit different temporal niches, thereby dominating in different seasons without outcompeting the other species (Figure 2e; Díaz et al., 1999). This may explain the positive effect sizes at high minimum temperatures and wide temperature ranges (Scheiner & Rey-Benayas, 1994;González-Moreno et al., 2015). Nevertheless, such a hypothesis should be investigated further by inspecting variation in species abundance in a community over time (Tilman, 1996;Pescador et al., 2015).
Further note that due to the disregard of intraspecific variation in this study, we potentially overestimate functional richness in environmentally constricted areas as adaptive traits may be more similar between species in these locations due to environmental filtering under the favorability hypothesis (Grant & Abbott, 1980;Swenson & Enquist, 2009;Swenson et al., 2012). Conversely, we may underestimate functional richness in other areas as adaptive traits may be more dissimilar between species under the principle of limiting similarity (Weiher & Keddy, 1995;Stubbs & Wilson, 2004;Mason & Wilson, 2006;Violle et al., 2011). Future field campaigns might invest in local trait measurements for all species to enable the inclusion of intraspecific trait variation (Albert et al., 2010;Niu et al., 2020).
Functional and taxonomic richness decreased toward higher minimum temperatures and narrow temperature ranges. Under the stress gradient hypothesis, the reduced diversity may indicate strong competition due to a benign, warm climate (Bertness & Callaway, 1994). However, temperature range is known to affect functional richness (González-Moreno et al., 2015) and might act as a stressor in these locations explaining the reduction in functional and taxonomic diversity under the favorability hypothesis. Since positive effect sizes were found in these environments, it could be a more stressful environment where facilitation possibly causes trait divergence.
This may follow the same explanation for the trait divergence found at extreme low minimum temperatures. In addition, facilitation was also observed in an alpine study where greater variation in temperature made for more stressful environments (Molenda et al., 2012).
Nevertheless, these opposite effects may suggest that the locations with high minimum temperatures and narrow temperature ranges might be restricted by additional stressors that were not included in this study, like summer temperatures, wind, solar radiation, or even different herbivory levels. Experiments are required to disentangle the causal mechanisms underlying the functional and taxonomic richness patterns. They may also determine if and which environmental factors determine the diversity in these European grasslands, but our results indicate the importance of climate seasonality.
Interpreting patterns in community diversity along climatic gradients is difficult, particularly because the variation in our data explained by the temperature and precipitation variables was relatively low (8-10%; Appendix S10). Results should be interpreted with caution and experiments are required to back up possible explanations of diversity patterns. The low explained variation indicates that other factors beside climate vary between locations and play a large role in determining the plant diversity of European grassland communities, such as variation in soil, herbivory, management, landscape history and various biogeographical influences (e.g., Willems, 1983;Bakker et al., 2006;Dainese et al., 2015). It should also be emphasized that functional diversity is always dependent on the traits that are included in the study, where results may change when different traits are selected (Villéger et al., 2008). This means that the conclusions regarding functional richness trends in this study depend on the specific set of traits we selected. Future research could benefit from the ever-growing availability of trait data (Kattge et al., 2020) and vegetation-plot data from other ecosystems (Bruelheide et al., 2019) to assess the generality of plant community diversity patterns along environmental gradients and to assess how community assembly may change in response to global warming (Mouillot et al., 2007;Mason et al., 2011;Catford et al., 2020).