Soil Abiotic Properties Shape Plant Functional Diversity Across Temperate Grassland Plant Communities

There is increasing awareness that plant community functional properties can be an important driver of ecosystem functioning. However, major knowledge gaps exist about how environmental factors, especially climate and soil abiotic properties, shape plant functional diversity at a regional scale. Furthermore, at those scales the relationships between plant functional and taxonomic diversity have rarely been considered. Here, we used a large database of plant species and functional trait data from 180 temperate grasslands across England, covering a broad range of grassland types, climatic conditions and management intensities, the last having a strong influence on multiple soil variables. Our specific aims were to: (1) identify the dominant environmental factors explaining variation in different facets of plant community functional properties, including community weighted means (CWMs) of functional traits and various multi-trait functional diversity indices; and (2) test whether the relationship between plant functional and taxonomic diversity is mediated by environmental factors at a regional scale. We found that soil abiotic properties (pH and nutrient stocks), but not climate, were the main environmental factors explaining grassland plant functional diversity at a regional scale, with a significant contribution of soil nutrient stoichiometry (N/P ratio). Two indices of plant community functional properties, namely CWMs of specific leaf area and relative growth rate, were explained by interactions between soil pH and N and mean annual precipitation, soil pH soil N and soil N/P ratio. These indices were also negatively related to taxonomic diversity under certain soil abiotic conditions, specifically high soil clay content, pH and N/P. Together, our results indicate that soil abiotic properties rather than climate factors shape plant functional diversity across temperate grassland plant communities at a regional scale. They also suggest that interactions between environmental factors play a significant role in shaping patterns of plant community functional properties. Our findings are of importance for the design and interpretation of future studies using trait and diversity measures as proxies of ecosystem services at regional scales.


ABSTRACT
There is increasing awareness that plant community functional properties can be an important driver of ecosystem functioning. However, major knowledge gaps exist about how environmental factors, especially climate and soil abiotic properties, shape plant functional diversity at a regional scale. Furthermore, at those scales the relationships between plant functional and taxonomic diversity have rarely been considered. Here, we used a large database of plant species and functional trait data from 180 temperate grasslands across England, covering a broad range of grassland types, climatic conditions and management intensities, the last having a strong influence on multiple soil variables. Our specific aims were to: (1) identify the domi-nant environmental factors explaining variation in different facets of plant community functional properties, including community weighted means (CWMs) of functional traits and various multi-trait functional diversity indices; and (2) test whether the relationship between plant functional and taxonomic diversity is mediated by environmental factors at a regional scale. We found that soil abiotic properties (pH and nutrient stocks), but not climate, were the main environmental factors explaining grassland plant functional diversity at a regional scale, with a significant contribution of soil nutrient stoichiometry (N/P ratio). Two indices of plant community functional properties, namely CWMs of specific leaf area and relative growth rate, were explained by interactions between soil pH and N and mean annual precipitation, soil pH soil N and soil N/P ratio. These indices were also negatively related to taxonomic diversity under certain soil abiotic conditions, specifically high soil clay content, pH and N/P. Together, our results indicate that soil abiotic properties rather than climate factors shape plant functional diversity across temperate grassland plant communities at a regional scale. They also suggest that interactions between environmental factors play a significant role in shaping patterns of plant community functional properties. Our findings are of importance for the design and interpretation of future studies using trait and diversity measures as proxies of ecosystem services at regional scales.

INTRODUCTION
The last 2 decades have witnessed a surge in interest in understanding how plant functional diversity regulates ecosystem functioning (Díaz and Cabido 2001;Moretti and others 2013;Cadotte 2017;Huang and others 2019). Functional diversity refers to those components of biodiversity that influence how an ecosystem operates (Tilman 2001). It is commonly measured using indices that account for the distribution of selected functional traits, which are morphological, physiological, phenological or behavioral features of an organism that can be measured at the individual level and that affect its fitness (Carmona and others 2016). For plants, some of the most commonly employed traits are those associated with the leaf economics spectrum (that is,, leaf nitrogen (N), leaf mass area (LMA), leaf lifespan and photosynthetic rates; Wright and others 2004). This describes a collection of leaf traits that reflect a trade-off between production of plant biomass and persistence of plant organs and individuals. This trade-off contrasts cost-effective short-lived leaves with rapid return on carbon (C) and nutrient investment, with more costly long-lived leaves with slow returns on investment (Reich 2014).
There are two main types of trait indices used to characterize community functional properties at a large spatial scale (Díaz and others 2007). The first type are functional dominance indices like community weighed means (CWMs), The first type are functional dominance indices like community weighed means (CWMs), which refer to the trait value most frequently expected in a community (Díaz and others 2007), and are used to test the mass ratio hypothesis, which states that ecosystem functioning is chiefly determined by trait values of the dominant contributors to plant biomass (Grime 1998;Díaz and others 2007). The second type of indices are functional diversity indices, that provide information about how trait values are distributed in the community (Díaz and others 2007). For instance, Rao's Q index is a multivariate form of variance (de Bello and others 2016) and functional divergence (FDiv) indicates if the trait values of the dominant species are centered or extreme within the functional trait range (Villé ger and others 2008). These trait distribution indices can be used to test for the niche complementarity effect, whereby combinations of species use resources complementarily, increasing resource use efficiency (Díaz and Cabido 2001). Clearly, both types of indices provide complementary information, since dominance and niche complementarity effects can both drive ecosystem functioning (Cadotte 2017; Zhou and others 2017).
The effects of many of these traits on ecosystem properties and processes, such as primary productivity, litter decomposition, soil nutrient cycling and soil carbon (C) or N storage, are well studied (Funk and others 2017 and references therein). However, although past studies have explored how leaf traits are related to climate and soil resource availability (Cunningham and others 1999;Poorter and De Jong 1999;Wright and others 2004;Maire and others 2015; Bruelheide and others 2018), it remains unresolved how different facets of functional diversity are shaped by environmental factors such as climate and soil conditions at large spatial scales, and which factors play a dominant role. In particular, the majority of studies exploring how environmental factors drive trait indices have focused on CWMs, omitting trait distribution indices (Dwyer and others 2015; Simpson and others 2016; Borgy and others 2017; Herná ndez-Vargas and others 2019; but see Barros and others 2018 and Busch and others 2018).
Another important knowledge gap is that the effects of some of the potential drivers of plant functional diversity indices, such as climate and soil properties, are partially unknown. This uncertainty exists for three reasons. First, most studies the drivers of plant functional diversity have studied a narrower range of values than the one that could be expected over larger environmental gradients. Hence, the potential effect of one factor could vary depending on the range of values considered for each variable (Maire and others 2015). Second, many studies have tested only a small set of potential environmental drivers or studied strong environmental gradients involving multiple highly correlated factors (Barros and others 2018; Herná ndez-Vargas and others 2019). This makes it difficult to distinguish between the effects of each considered factor and implies that interactions be-Soil Abiotic Properties Shape Plant Functional Diversity Across Temperate Grassland tween potential drivers of functional diversity indices are not frequently tested. Third, some potential drivers of plant functional diversity, such as soil nutrient stoichiometry, are rarely considered despite being an important driver of plant leafeconomic traits (Fujita and others 2014); separating nutrient imbalance limitation effects is crucial for understanding ecological processes (Cardinale and others 2009).
Relationships between plant community functional properties and taxonomic diversity are also of importance for biodiversity conservation since the dominance and/or the variability of some plant traits in some ecosystems could be more related to species diversity and loss than others (Wardle 2016). Furthermore, some indices of plant community functional properties are closely related to taxonomic diversity indices and affected by similar drivers, while others may be not be closely coupled (Zhou and others 2016). Strong coupling is expected between Rao's Q and species richness, since the more plant species there are in the community, the more diverse the expected trait pool (Petchey and Gaston 2002;De Bello and others 2013). However, relationships are more heterogeneous when other indices of plant community functional properties are considered, especially CWMs (Zhou and others 2016; Busch and others 2018). Additionally, functional-taxonomic diversity relationships may be mediated by environmental drivers (Cadotte and others 2011;Morelli and others 2018). However, to our knowledge, this has not been tested at large spatial scales, such as across plant communities within regions, nor for a wide range of plant trait indices.
Here, we identified the dominant potential abiotic drivers shaping different facets of plant functional diversity and examined how functional diversity-taxonomic diversity relationships vary with environmental drivers across a broad range of temperate grassland plant communities at a regional scale. We addressed two specific questions: (1) what are the dominant environmental factors shaping CWM traits and multi-trait functional diversity indices across temperate grassland plant communities? and (2) is the relationship between taxonomic diversity and functional diversity indices mediated by environmental factors? This was tested at a regional scale using a large dataset of plant species and functional trait data and associated information on a range of soil abiotic properties and climatic variables. This dataset comprised 180 grasslands across England, subject to different levels of historical management intensity (Table S1) and covering a broad range of grassland habitats and climatic and soil conditions. Sampling campaigns were carried out during the summer of 2005, which presented few temperature and precipitation anomalies (Met Office 2022). These data have been used previously to examine relationships between plant CWM traits and surface soil properties, including microbial community composition and diversity (de Vries and others 2012; Delgado-Baquerizo and others 2018) and soil C stocks at a regional scale (Manning and others 2015). We advance on this work by testing the relative importance of environmental factors, including soil abiotic properties and climate, on different facets of plant functional diversity and plant functional-taxonomic diversity relationships at a regional scale.

METHODS
A wide range of grasslands from twelve regions of England ( Figure S1) including acid, calcicolous, mesotrophic and wet grasslands (Rodwell 1992), and covering a broad range of soil and climate conditions, were sampled during June and July 2005. Within each region, there were five sites, each containing three fields with three different levels of historical management intensity: unimproved (U) grassland, semi-improved (SI) or improved (I) grassland, resulting in a total of 180 fields. Unimproved grasslands are defined as having no history of fertilizer application and are typically species-rich and of high conservation value, whereas semi-improved and improved grasslands have a history (> 10 years) of receiving 25-50 kg N ha -1 year -1 and > 100 kg N ha -1 year -1 , in the form of inorganic fertilizer, respectively. Additionally, improved grasslands are typically subject to higher grazing pressures and/or more frequent cutting (Table S1; Ward and others 2016). The three fields at each location were selected to be of a similar soil type and topography, so it was assumed that the only major difference between fields was historic management. According to the UK National Vegetation classification (Rodwell 1992), the semi-improved and improved grasslands were mainly the species-poor Lolium perenne-Cynosurus cristatus (MG6) and Lolium perenne (MG7) communities, respectively. Conversely, there are many different, but typically species rich, plant communities of unimproved grasslands (Ward and others 2016). Within each field, a 25 9 25 m plot of homogenous vegetation was determined, estimating the percentage cover of each plant species in three 1 m 2 quadrats placed at random. In each quadrat, five 2-cm-diameter 7-cm-deep soil cores were taken at random and pooled to produce a composite sample for chemical analysis. Soils were sampled to 7-cm depth following the UK's Department of Environment, Food and Rural Affairs (DEFRA) recommended sampling depth for assessment of soil abiotic properties in permanent grassland (DEFRA 2010).
Gravimetric moisture content was determined, and air-dried soil was analyzed for total N, phosphorus (P) and pH, using standard methodology. Total soil N and P were measured on an Elementar Vario EL elemental analyzer (Hanau, Germany). A standardized wet sieving method was used (De Deyn and others 2011; Manning and others 2015) to separate the soil particle size fractions. Two sieves were used (250-and 50-lm meshes). Soil texture classification was based on the proportions of sand, silt and clay using the triangular classification of the Soil Survey of England and Wales RDS Technical Advice Note 52 (Anonymous 2006). Bulk density was calculated after removal of all stones and roots > 3 mm diameter and this value was used to calculate soil measures on a per unit area basis. More information about the analyses of soil used here can be found in de Vries and others (2012) and Manning and others (2015). Hereon, we describe all the variables used in this study (Table 1). An overview of the range and average values of soil properties and plant diversity indices is given in Table S2.

Climate Data
Each field was assigned to a 5 km 9 5 km grid containing long-term gridded climate data obtained from Met Office UKCP09 databases (Jenkins and others 2009). Mean annual temperature (MAT) and mean annual precipitation (MAP) were calculated from monthly temperature and rainfall averages from 1981 to 2006, ranging between 6.5-10.9°C and 596-3191 mm, respectively.

Plant Taxonomic Diversity and Trait Functional Properties Data
We estimated plant taxonomic diversity on each site employing plant species richness, defined as the number of species recorded per site, and Camargo evenness index (Camargo 1993), calculated with the following formula: where p i1 and p i2 are the relative abundance of the i1th and the i2th species, respectively. We included an evenness index since it complements information from richness estimates (Stirling and Wilsey 2001; Wilsey and others 2005; Hillebrand and others 2008), and we selected this index as it is not a function of species richness (Tuomisto 2012). In addition, this variable was less redundant in rela-  tion to species richness according to the variance inflation factor than other evenness estimations (Miles 2005).
For estimating different indices of plant community functional properties, we selected four commonly used plant traits: dry matter content (LDMC), specific leaf area (SLA), relative growth rate (RGR) and leaf nitrogen content (leaf N). For consistency with previous works (de Vries and others 2012; Manning and others 2015), LDMC values of each species were obtained from Grime and others (2007), whereas the rest trait values were obtained from the TRY database in the context of earlier analyses (Kattge and others 2011; de Vries and others 2012; Manning and others 2015). For some species, TRY records were dominated by a single research project or contributor. To prevent biasing species means toward those of their specific methodology or study system, trait values were therefore averaged within each contributing author and then across authors to obtain species level trait values. The main data sources contributing to the TRY dataset we used for obtaining these trait values were others (2003, 2004), Garnier and others (2007), Kleyer and others (2008), Shipley (1989Shipley ( , 2002, Shipley and Vu (2002) and Shipley and Parent (1991). We used trait data from databases instead of direct trait measures since with this approach we avoid the introduction of noise from factors like individual and phenological variation, intraspecific trait variation is not accounted for but is typically small relative to total trait variation when using averaged weighted values, as done in this study (Kichenin and others 2013). Thus, the average trait values we used for each species represents an average value across a range of environmental conditions. These values were considered meaningful because interspecific variation is much greater than intraspecific variation for the traits selected here (Kattge and others 2011).
We calculated community weighted means (CWMs) of each trait as: Being p i the relative contribution of species i to the cover of n which pi is the relative contribution of species i to the cover of the whole community, n is the number of most abundant species, and trait i is the trait value of species i (Garnier and others 2004).
We used the R package FD (Laliberté and others 2015) to calculate two multi-trait distribution indices: Quadratic entropy (Rao 1982;Botta-Duká t 2005), and Functional divergence (FDiv; Villé ger and others 2008). Rao's Q expresses functional diversity as the average of the species differences considering a measurement of pairwise differences between species and relative frequency data. FDiv quantifies how the trait values are spread along the range of a trait space.

Statistical Analyses
All statistical analyses were performed with the software R ver. 3.4.3 (R Core Team 2017), at 95% significance level. All models fitted were linear mixed models with a site random effect to consider the clustering of fields in triplets. We also carried out mixed-effects AVOVA models to assess how management intensity shaped soil variables and diversity indices (Table S2).
We fitted models for explaining indices of plant community functional properties with abiotic factors and plant taxonomic diversity factors. These models aim to characterize the effects of the abiotic factors and their relationship with taxonomic diversity indices, including the detection of interaction effects. Additionally, we generated models for explaining taxonomic diversity variables to identify the abiotic factors affecting both plant community functional properties and taxonomic diversity indices. The results of the models for taxonomical diversity indices can be found in Appendix.
Before starting, we checked for curved relationships by exploratory methods like the scatterplot function of the car package (Fox and others 2018) and by building single predictor models, considering the squared term of the variable or applying a log transformation. If including squared or logtransformed terms improved the model fit, we considered those terms in the modeling procedure. Hence, our modeling procedure took into account curved effects of predictors on the response variables, (de Vries and others 2012).
We used a modified version of the 'hierarchical' modeling procedure employed in previous studies based on this database (de Vries and others 2012; Manning and others 2015) and following Díaz and others (2007). A variation of this method can be found in Rodríguez and others (2020). We grouped functional diversity drivers in four groups: climate, soil type, soil stocks and taxonomic diversity variables (Table S2). Climate variables were MAP and MAT. Soil type variables included soil physicochemical properties, namely soil pH, moisture content and soil texture and soil stocks, which included soil nutrient and C stocks (total N, P and C, and C/N and N/P ratios). Taxonomic diversity variables were plant species richness and evenness. We added driver groups according to a fixed sequential order, selecting model terms automatically by a genetic algorithm included in the R package glmulti (Calcagno 2015). With this procedure, we ensured that the variables that are the ultimate cause of variation of each response variable were added before other predictors affected by these variables. We did not include management intensity as a potential predictor, since it is a classification based on inorganic nitrogen supply and its inclusion could mask the relevance of soil nutrients and C stocks as drivers of plant functional diversity, and also increase the complexity of the results, hampering their interpretation. However, the differences in measured soil variables along the management intensity gradient are given in Table S2. We also created models including management intensity between soil type and soil stocks variables in the hierarchical modeling procedure (de Vries and others 2012) for all functional diversity and present these in Appendix. The criterion for model selection was the corrected Akaike information criterion (AICc), which is appropriate when n/k is less than 40, being n the sample size and k the number of parameters in the most complex model (Symonds and Moussalli 2011). The detailed methodology of each of these steps is fully described in Appendix.
Once we had the final model, we assessed the significance of each term by removing it and performing a LRT. For estimating the significance of the main effects, we also removed the interaction terms in which they were involved to avoid transferring the effects of the main terms to the interaction terms. We estimated the marginal variance explained by the fixed effects by the method of Nakagawa and Schielzeth (Nakagawa and Schielzeth 2013), using the r.squaredGLMM function from the R package MuMIn (Barton 2018). The importance of each potential driver groups (climate, soil type, soil stocks, and plant taxonomic diversity) was estimated as the difference between the total marginal variance explained by the model and the marginal variance explained by a model where the variables belonging to the corresponding group were suppressed.
To look for significant differences between variable group importance values (climate, soil type, soil stocks and plant taxonomic diversity), we estimated 95% bootstrapped confidence intervals using the package boot (Canty and Ripley 2017). We used the bias-corrected accelerated method (BCa) as several bootstrapped distributions were biased (Hesterberg and others 2005). We also determined the importance of each term in the model by observing AICc change (D) on deletion.
As a measure of trait variance, a larger Rao's Q could just be caused by a larger species pool (Barros and others 2018). We refitted Rao's Q model substituting it by a standardized effect size (SES) measure calculated with the package picante (Kembel and others 2010) following de Bello and others (2016). With this, we tested if the effects driving Rao's Q remained when removing the species richness effect (Barros and others 2018).

RESULTS
Historical grassland management was a strong driver of soil abiotic properties (Table S2), with a trend of decreasing soil moisture and C/N and N/P ratios, and increasing total N and total P along a gradient of increasing management intensity. Soil pH did not show significant differences across intensity levels. Additionally, most plant diversity indices varied with management intensity: species richness, LDMC, Rao's Q and FDiv decreased across the gradient of increasing management intensity, whereas LeafN, SLA and RGR increased. Average values on each management type of these variables are shown in Table S2.

CWM Plant Traits
The fixed effects of the CWM plant trait models explained between 46 and 58% of the variance and their residuals did not show important deviations from the normal distribution (Appendix 1). The models obtained for CWM LDMC and leaf N were simpler than those for RGR and SLA, with fewer drivers and interaction terms (Tables 2, 3, S3-S6).
Soil type and nutrient stocks were the two variable types explaining LDMC ( Figure S2), and with similar importance according to the bootstrapped confidence intervals ( Figure S2). The most important explanatory variables in the LDMC model were soil pH, N/P ratio and total N (27.9, 45.6 and 11.3 DAICc values; Tables 2 & S3). LDMC responses to soil pH and total N had inverse curve shapes with a minimum value of 6 ( Figure 1B) and 0.8 kg m -2 ( Figure 1C), respectively. Conversely, LDMC peaked at an intermediate N/P ratio value of 15 ( Figure 1D). Additionally, LDMC decreased with increasing soil clay content.

Soil Abiotic Properties Shape Plant Functional Diversity Across Temperate Grassland
Leaf N wan explained by soil type, soil nutrients and plant taxonomic diversity with no significant differences between groups according to the bootstrapped confidence intervals ( Figure S2 Figure 3B). SLA values increased with total soil N under low pH conditions ( Figure 3B). Additionally, SLA decreased with increasing soil N/P values ( Figure 3A).
Finally, RGR was affected by climate, soil type, soil nutrients and plant diversity variables. Bootstrapped confidence intervals identified significant differences between plant diversity and climate, accounting for soil type and nutrients at intermediate values ( Figure S2). Species richness was the most important variable for RGR (DAICc = 54.5), followed by soil pH and the N/P ratio (31.7 and 25.4 DAICc values, respectively; Tables 2 & S6, Figure 4). Soil pH produced a maximum peak on RGR under low MAP and plant species richness values, with a decreasing effect with increases in these variables ( Figure 4A, D).

Multi-trait Functional Diversity Indices
The fixed effects of the Rao's Q and FDiv models explained between 57 and 38% of the variance, respectively, and their residuals did not show important deviations from the normal distribution (Appendix III). Unlike plant diversity and CWM plant trait models, these two models included total soil C as a significant driver (Tables 2 and S7-S8). Rao's Q was driven mainly by plant diversity (Figure S3), which explained 24% of the variance, while soil type and stocks played minor roles. Soil pH and total P were the most important soil drivers of Rao's Q (DAICc = 21.13 and DAICc = 29.88, respectively). Rao's Q reached its minimum mean value at the intermediate pH value of 6 and decreased when soil total P increased ( Figure 5A and C). When comparing Rao's Q model with the SES model ( Figure S4), only the relationship with taxonomic diversity was changed, indicating that soil abiotic properties were only acting on Rao's Q. Additionally, Rao's Q increased with increasing soil C values ( Figure 5B) and reached its maximum mean value at the intermediate evenness value of 0.5 ( Figure 5E).
FDiv was mostly related to soil nutrient stock variables, which explained 34% of its variance ( Figure S3). Soil N/P ratio was the most relevant variable for FDiv (DAICc = 4)0.88; Table 3), with the N/P effect on FDiv being positive at low evenness and negative at high evenness values (Figure 6C). Furthermore, Fdiv increased with soil C ( Figure 6A).  Table 2.

Functional-Taxonomic Diversity Relationship
CWM LDMC was the only functional trait measure that was not related to taxonomic diversity. Plant species richness was negatively related to leaf N ( Figure 2D), but interacted with several abiotic variables in the SLA and RGR models. Species richness was associated with lower SLA in soils with high clay values ( Figure 3C) and low RGR values with high pH and N/P levels, but this effect disappeared at low soil pH and N/P ( Figure 4D and E).
Plant species richness was by far the most important variable for Rao's Q (DAICc = 43.99; Table 3), exerting a positive effect ( Figure 5D). Rao's Q also was related to plant species evenness (DAICc = 25.9; Table 3), peaking at the intermediate value of 0.5 ( Figure 5E). FDiv was related to plant species evenness (DAICc = 13.83; Table 3) in interaction with the soil N/P ratio: FDiv was positively related to plant species evenness at low N/P, but negatively at high N/P ( Figure 6E).

DISCUSSION
In this study, we used a large database compiling plant, soil and climatic information from 180 grasslands across England, spanning all major grassland types and different levels of management intensity, to find the dominant environmental factors shaping plant functional diversity, and to characterize the relationship between plant functional and taxonomical diversity at a regional scale. Figure 2. The relationships between leaf N and climate, soil type, soil nutrient and plant diversity predictors. In ''a-c'' the solid line represents the relationship between leaf N (y-axis) and the corresponding predictor (x-axis). In ''b'' the relationship shown is leaf N explained by S and E. All the predicted values are shown keeping all other predictors at their average value. For statistics, Table 3.
We obtained a total of six models that explained a noticeable portion of the variance in plant functional properties across English temperate grasslands. The main highlights of these models were that first, soil abiotic factors, including nutrient stoichiometry (that is, N/P ratio), were the main factors explaining plant functional diversity across grasslands at a regional scale, rather than climatic factors; and second, the relationship between plant functional and taxonomical diversity is mediated by soil abiotic factors at a regional scale. As expected, variation in soil abiotic properties across temperate grasslands was influenced by historical grassland management, with most soil variables varying along the gradient of management intensity (Table S2). As such, our findings underscore the important role of grassland management in shaping grassland plant functional diversity via its impact on abiotic properties related to the availability of nutrients in soil.

Drivers of Plant Trait Functional Diversity
Overall, soil pH, total N stock and nutrient stoichiometry (that is, soil N/P ratio) were the most important variables explaining plant functional diversity indices. Soil pH was the most important soil variable in all CWM models, showing a humpshaped relationship with leaf N, RGR and SLA models, and the inverse pattern in the LDMC model. This relationship is likely explained by intermediate pH values maximizing plant nutrient availability, because low or high pH values affect soil fertility by reducing the availability of different Figure 3. The relationships between SLA and climate, soil type, soil nutrient and plant diversity predictors. In ''a'' the solid line represents the relationship between SLA (y-axis) and N/P (x-axis). In ''b-c'' the relationship shown is SLA explained by a pair of interacting predictors. All the predicted values are showed keeping all other predictors at their average value. For statistics, Tables 3 and S3. macro-and micronutrients (McGrath and others 2014). For instance, soil P is strongly adsorbed by mineral surfaces at low pH values, but can precipitate as Mg and Ca phosphates at high pH values (McGrath and others 2014).
It is generally expected that plant communities dominated by species with resource-exploitative traits (that is, high leaf N, SLA and RGR) tend to occur on soils with an abundance of resources, whereas communities dominated by species with conservative traits (that is, high LDMC) predominate on soils with low available nutrients (Reich 2014; Allan and others 2015; Boeddinghaus and others 2019). Accordingly, we found that intermediate soil pH values were also associated with lower plant Rao's Q ( Figure 5A), suggesting that  Tables 3 and S4. niche space is reduced through the displacement of nutrient-conservative poor light competitor species (Hautier and others 2009). In addition, soil pH appeared to alter the MAP effect on RGR, which was greater at higher MAP values under low pH conditions, while at medium-high pH values the relationship between MAP and RGR was negative.
This effect could be related to the peak in soil bacterial biomass found by de Vries and others (2012) under high MAP and low pH values, which could promote a faster nutrient turnover and supply (Grigulis and others 2013).
Considering the important role of soil pH in determining nutrient availability, the role of sev- Figure 5. The relationships between Rao's Q and soil type, soil nutrient and plant diversity predictors. The solid lines represent the relationship between Q (y-axis) and the corresponding predictor (x-axis). All the predicted values are showed keeping all other predictors at their average value. For statistics, Tables 3 and S5. eral soil nutrient variables played in driving plant functional diversity was not unexpected (Figure S2). However, the soil N/P ratio, which decreased with increasing grassland management intensity, was one of the most important soil variables in several models, supporting the hypothesis that not only nutrient availability, but also nutrient stoichiometry, shapes plant functional diversity (Fujita and others 2014). Low N/P environments, typical of more intensively managed grasslands that receive high rates of N fertilizer relative to P (Tables 2 and S1), favor nutrient exploitative plant species, whereas high N/P environments favor nutrient conservative plant species (Fujita and others 2014). For instance, both low total N and high N/P, were associated with high LDMC and low leaf N, in accordance with the leaf economics hypothesis (Wright and others 2004). CWM SLA was also higher under low N/P conditions, whereas total soil N effects on SLA were only positive in acidic soils. Provided that available N losses via ammonia volatilization mostly occur under pH values > 8 (above those considered in our study (Freney and others 1983;McGrath and others 2014), the lack of relationship between soil N content and SLA at high pH values might be due to precipitation of other nutrients, such as phosphorus, under alkaline conditions (McGrath and others 2014).
The CWM RGR model also highlights the importance of soil nutrient stoichiometry as a driver of plant functional composition, since the soil N/P ratio had several interactive heffects (Figure 4). Fast growing plants are commonly promoted where there are no nutrient and water limitations (Aerts and Chapin 1999; Fujita and others 2014). However, our results point a decreasing importance of soil N/P for RGR in higher Figure 6. The relationships between FDiv and soil nutrient and plant diversity predictors. In ''a-b'' the solid line represents the relationship between FDiv (y-axis) and the corresponding predictor (x-axis). In ''c'' the relationship shown is FDiv explained by a pair of interacting predictors: N/P and E. All the predicted values are showed keeping all other predictors at their average value. For statistics, Tables 3 and S6. rainfall sites, and an increasing effect of N/P at high soil total N values. RGR is a good estimator of plant productivity (Vile and others 2006), and consequently, the complexity of the RGR model suggests that the productivity of temperate grassland communities is highly variable and depends on multiple environmental factors.
As was the case for pH, high-soil-nutrient conditions (soil P) and nutrient imbalance (N/P) were negatively related to functional diversity in multitrait distribution indices, thus supporting the hypothesis that niche space is reduced in nutrientrich soils (Harpole and Tilman 2007;Hautier and others 2009). Interestingly, soil C content was selected only in these last two models. Soil C constitutes an indicator of soil organic matter as its main component (Ondrasek and others 2019). Hence, the positive relationship of soil C content with Rao's Q and FDiv is potentially indicative of the multiple effects of soil organic matter on soil biochemistry (Ondrasek and others 2019), which may increase plant niche availability without affecting plant functional dominance.
Climatic variables were rarely selected in our models, despite the study sites covering a wide range of MAT and MAP (6.5-10.9°C and 596-3191 mm, respectively). Only RGR was related to MAP, decreasing at high MAP values under high pH and N/P ratio values. This result is not unexpected and can be understood via several lines of evidence. First, many leaf strategies are viable in most climate conditions (Reich 2014), which often leads to weak correlations between plant traits and climate (Wright and others 2004;Reich 2014). Second, plant trait responses to climatic variables are heterogeneous (Griffin-Nolan and others 2018) and depend on the climate regimes and plant species and/or traits considered (Moles and others 2014; Borgy and others 2017; Schellenberger Costa and others 2017). Third, broad scale climate factors, taken from climate models, can be of modest importance in models in comparison with plotmeasured soil variables (Bruelheide and others 2018). Furthermore, there are at least two possible sources of underestimation in the importance of climatic drivers on functional diversity in this study. First, climate can condition plant functional diversity exclusively through the remaining drivers, and hence its effects could be masked by other predictors (Maire and others 2015). Second, simple temperature and rainfall variables like MAP and MAT could be poorly related to plant trait indices compared with indices of seasonal or growing season length (Borgy and others 2017).

The Plant Functional-Taxonomic Diversity Relationship
Our models showed that plant functional-taxonomic diversity relationships can vary depending on traits and the indices considered, sometimes involving interactions with environmental drivers. Despite controlling for climate and soil factors in our models, plant species richness was included in all CWM plant trait models except the model for LDMC, being the most important single variable for predicting leaf N, RGR and SLA. The absence of this predictor in the LDMC model could be due to the fact that high or low LDMC values are actually strategies related to stability of plant populations over time (Má jekova and others 2014) rather than dominance processes linked to local extinctions.
The inclusion of plant species richness in the rest of our CWM models provided us information about which traits are selected under species rich and poor conditions independently of climate and soil factors. The negative relationships between species richness and CWM leaf N, SLA and RGR are likely explained by slow growing nutrient-conservative species being outcompeted by fast growing resource-exploitative species (Saar and others 2012), which is commonly reported in grasslands subject to nutrient enrichment (Stevens and others 2004;Borer and others 2014). CWM SLA remained stable across the species richness gradient in low clay content soils. This pattern could be due to the high homogeneity, in terms of plant ecological niches, of soils with low clay content (Laureto and Cianciaruso 2015).
The relationship between RGR and plant species richness was only negative under high-soil-pH and N/P conditions. This interaction of species richness with N/P determining RGR supports the hypothesis that nutrient-rich conditions enhance competition for light (Hautier and others 2009). In grasslands already dominated by a few high RGR species, there was no N/P effect detected, but in species-rich grasslands low N/P would release light competition, excluding low RGR species. Provided that RGR only shifted to low values in species-rich grasslands on alkaline soils, the model suggests that at high pH values there would be a decrease of nutrient availability promoting slow growing nutrient conservative plant species (McGrath and others 2014; Reich 2014).
Rao's Q was related to both plant species richness and evenness. The observed linear relationship between Rao's Q and species richness detected was unexpected: the relationship is typically positive at low species richness values and neutral at high Soil Abiotic Properties Shape Plant Functional Diversity Across Temperate Grassland species richness values (Petchey and Gaston 2006;de Bello and others 2016). Pavoine and others (2013) suggested that strong correlations between functional and taxonomic diversity may be related to the existence of a stress gradient, which suggests that there may be important unmeasured factors. Additionally, the positive relationship between species richness and Rao's Q also indicates that each species could contribute to ecosystem functioning via its particular set of functional traits (Naeem and Wright 2003). On the other hand, Rao's Q peaked at intermediate evenness levels, so dominant species could have a role in trait assemblage (Hillebrand and others 2008).
The relationship between evenness and FDiv was modulated by soil N/P ratio, which, as mentioned earlier, decreased with increasing management intensity. This might be explained by high evenness being associated with high light use efficiency and complementarity of plant communities (Wang and Yu 2018), and/or increasing functional divergence at low N/P conditions, when competition for light intensifies (Hautier and others 2009). In low-N/P conditions, plant divergence was negatively associated with evenness, probably because trait strategies of dominant and subordinate species diverge (Grime 2006).

CONCLUSIONS
Our findings indicate that soil abiotic properties, which are influenced by management intensity, are the main drivers of functional diversity across temperate grassland plant communities at a regional scale. Most of these effects could be attributed to the fast-slow plant economics spectrum theory, with more nutrient acquisitive traits being associated with higher soil nutrient availability. For instance, intermediate soil pH values, which are typically associated with high nutrient availability, were consistently associated with nutrient exploitative plant traits across the considered plant diversity indices. Moreover, our models suggest that not only soil nutrient stocks, but also soil nutrient stoichiometry (N/P), acts an important driver of plant functional diversity patterns across temperate grasslands. In contrast, climate factors were of less importance in explaining plant functional diversity indices. Additionally, our results point to the importance of interactions between environmental variables in explaining plant functional diversity patterns, which are not usually considered in studies of the drivers of plant functional diversity. Finally, in SLA, RGR and FDiv models the relationship between plant functional and taxonomic diversity was modulated by abiotic soil properties. This indicates that the plant functional-taxonomic diversity relationship not only varies depending on the traits and indices considered, but also on environmental conditions. This finding should be contemplated in further studies addressing the relationship between plant traits and species loss. Although our results are based on plant trait values from databases (Kattge and others 2011), they provide a strong basis for future studies using functional and taxonomic diversity measures as proxies of ecosystem services at large spatial scales.