Environment and space drive the community assembly of Atlantic European grasslands: Insights from multiple facets

Ecological communities are assembled by regional and local processes. These processes select species through their traits, which are tied to species' evolutionary history. A multifaceted approach, encompassing taxonomic, functional and phylogenetic diversity can thus help us to better understand community assembly. We asked what the relative importance of geography, climate and soil parameters is in driving taxonomic, functional and phylogenetic diversities at local (alpha) and larger scales (beta diversity).


| INTRODUC TI ON
Understanding how plant species are assembled into communities is a longstanding goal in ecology. According to community assembly theory, plant species are assembled from region to local sites due to a combination of stochastic and deterministic processes, mostly related to geography, environment and/or biotic interactions within sites (Götzenberger et al., 2012;Kraft et al., 2015;Swenson et al., 2011). For example, by increasing geographical distance among sites, species diversity and composition might change because some species are not able to colonize more distant sites (i.e. dispersal limitation) or are not part of the available pool of species for a given locality (Chase & Myers, 2011). Likewise, geographically more distant communities might also be dissimilar in terms of environmental conditions (i.e. abiotic filter) and biotic interactions (e.g. biotic filters: competition, predation, herbivory etc.), filtering out species not able to cope with prevailing environmental conditions or tolerate certain biotic interactions at the local site (Götzenberger et al., 2012;HilleRisLambers et al., 2012;Kraft et al., 2015).
Disentangling the relative importance of these filters is challenging because they often operate simultaneously and their importance may change depending on the spatial scale and environmental gradients (e.g. climate, soil and water; Cornwell & Ackerly, 2009;Kraft et al., 2015). For example, climate is expected to be a large-scale environmental filter species have to overcome, whereas soil conditions are expected to be more important driving species diversity and composition at local scales (de Bello et al., 2013). In addition to this scale issue, the filtering process can affect multiple biodiversity facets (i.e. taxonomic diversity [TD], functional diversity [FD] and phylogenetic diversity [PD]), both within (alpha diversity) and among (beta diversity) sites. Therefore, considering different biodiversity facets, linked to spatial and environmental gradients, might provide complementary and detailed information on the plant community assembly processes.
Alpha and beta TD provide important insights regarding how the filtering process affects the number of species and their turnover (Cao et al., 2019). However, the potential of TD to explain community assembly and turnover is rather limited because it provides information on where species are, but not why. Instead, FD and PD can be used to understand why and how species are filtered within local sites. Functional traits are linked to the species' ability to cope with environmental conditions and competition, whereas PD can be used as a proxy and complement to FD, since closer lineages tend to possess similar traits and measuring all important traits in a community is not feasible (de Bello et al., 2017). FD can be calculated based either on multiple or single traits, enabling testing both the overall and individual trait contribution to community assembly patterns (Lepš et al., 2006). Thereby, considering simultaneously alpha and beta FD and PD, in addition to TD, allows to understand community assembly by showing which plants are selected by given environmental conditions based on their traits and shared evolutionary history (Pavoine & Bonsall, 2011;Villéger et al., 2013). Likewise, the effect of environmental and biotic filtering can be decoupled between alpha and beta diversities. For example, it is possible that a study region exhibits a high beta TD (different species identities) but low beta functional or phylogenetic diversities (similar functions and lineages), being an indication of strong environmental filtering (de Bello et al., 2009;de Paula et al., 2020). Although considering TD, FD and PD, at both alpha and beta levels, is key to disentangle the effect of multiple filters at different scales, few studies have used such a multifaceted approach to understand the community assembly of plants (but see, Dainese et al., 2015;de Paula et al., 2020;Swenson et al., 2012).
Here, we focused on the effects of geography, climate and soil in shaping Violion caninae grassland communities through multiple biodiversity facets, on both alpha and beta levels. European grasslands are an excellent model system to assess community assembly processes through a multifaceted approach due to their large geographical extent and extensive data sampling efforts, allowing to examine community assembly at both smaller and larger scales.
We used available data presenting 153 acidic grassland sites across Atlantic Europe, spanning from the South-West of France to the North of Norway (Stevens, Duprè, Dorland, et al., 2011). Combining these vegetation data with recent functional trait, phylogenetic and environmental data made the integration of multiple diversity facets (TD, FD and PD) at both alpha and beta levels feasible (Durka & Michalski, 2012;Kattge et al., 2020;Stevens, Duprè, Dorland, et al., 2011), providing a unique setup to test community assembly. Therefore, we expect that (1) taxonomic and phylogenetic beta diversity will be high due to the large geographical gradient examined, but functional beta diversity will be low because climatic conditions of the Atlantic region might act as a large-scale filter, sorting species with similar characteristics. Despite this, we expect that a considerable amount of variation in the diversity facets will occur.
In particular, (2) variables spanning larger gradients (i.e. geography and climate) should be key for beta diversity, while variables that are expected to be more heterogeneous at smaller geographical scales (i.e. soil and management parameters), might be relatively more important for alpha diversity. Finally, (3) as soils vary in nitrogen and environmental gradients is key to understand the relative importance of multi-scale processes driving the community assembly of European grasslands.

K E Y W O R D S
alpha diversity, assembly rules, beta diversity, biodiversity facets, biotic interactions, dispersal limitation, environmental filtering, equalizing fitness, functional redundancy phosphorus content across the study region (Helsen et al., 2014;Stevens, Duprè, Dorland, et al., 2011), plants will compete more strongly for these resources in nutrient rich sites, leading to the exclusion of less competitive species from the community. Therefore, at the alpha diversity level, nutrient rich sites might lead to lower TD, FD and PD, under the assumption that more competitive species share similar functional characteristics and/or are closely related (i.e. equalizing fitness processes; Chesson, 2000;Mayfield & Levine, 2010). Since the FD facet can be assessed through both multiple and single traits, we considered how both the overall alpha FD (all traits) and FD of single traits (seed mass, height and specific leaf area) vary along gradients.

| Vegetation data
A dataset of 153 acid grasslands belonging to the V. caninae association (Schwickerath, 1944) was used as provided by Stevens, Duprè, Dorland, et al. (2011), covering the Atlantic biogeographic zone of Europe and including nine countries ( Figure 1). All plots were surveyed between May and September in 2002 and 2007 and consisted of unfertilized grasslands that were managed by grazing or cutting. To standardize the data sampling across the geographical gradient, Stevens, Duprè, Dorland, et al. (2011) used a list of 14 indicator species that had to be found on a site. At each site, five randomly located 2 × 2 m quadrats were surveyed within a 1 ha area, estimating coverage of all vascular plants and bryophytes species using the Domin scale (values between 1 and 10). Because soil, management and climatic variables were only available at the site level, for this study, we pooled together the five quadrats at each site and averaged species abundances to obtain one aggregated sample for each site.

| Environmental data
At each site, 12 soil variables were sampled: pH and concentration of metals, nitrate, ammonium, total C and N, and Olsen extractable phosphorus (Stevens, Duprè, Dorland, et al., 2011; Appendix S1, Table S1). Measures of soil pH were available for soil samples from two different depths, but we only used values for the topsoil layer as they were highly correlated with those of the lower soil layer (r = 0.86, p < 0.001). Moreover, every site was described by inclination, vegetation height and two management variables: management type (grazed or mown) and grazing intensity (ordinal scale from 0 to 3; see Stevens, Duprè, Dorland, et al., 2011). Climatic variables were obtained from the CHELSA database (Karger et al., 2017). In total, six bioclimatic factors were considered in our study that are all averages of records over the period from 1970 to 2000, thus describing the prevailing recent climate: annual mean precipitation, isothermality, temperature annual range, annual precipitation, precipitation of driest month and precipitation seasonality (Appendix S1, Table S1).
To summarize climatic, soil and additional variables, we performed centred principal component analysis (PCA), separately for climatic, and the soil and management variables. Site scores along the first two PCA axes were used as variables reflecting either climatic conditions (Climate PCA) or soil and management properties (Soil PCA) of the sites.
Climate PCA axis 1 (explaining 48.5% of variance in climatic variables) was mainly represented by a precipitation gradient and precipitation seasonality, with sites with high precipitation also presenting high precipitation seasonality (Appendix S2, Figure S2.1). Climate PCA axis 2 (20.7% of variance explained) corresponded to temperature (mean annual temperature) and isothermality, with warmer sites also being more isothermal (i.e. diurnal temperature differences more similar to annual differences; Appendix S2, Figure S2.1). Soil PCA axis 1 (28.1% of variance explained) reflected a gradient of pH, general nutrient level (including ammonium) and land use (from ungrazed/mown to intensely grazed), while soil PCA axis 2 (10.4% of variance explained) mainly represented a nitrogen and nitrate gradient (Appendix S2, Figure S2.2). Additional PCA axes were not considered since their inclusion did not lead to better fitting statistical models (see Section 2.6).

| Geographical structure of study sites
To represent the geographical structure of the sampled sites and to account for spatial autocorrelation, geographical variables were transformed into site scores along the axes of the Moran's Eigenvector Map (MEM) of the geographical coordinates (Dray et al., 2006). MEM maps create spatial vectors by computing eigenvectors of a connectivity matrix based on the geographical position of the sampling sites (Griffith & Peres-Neto, 2006). The first spatial vectors represent a larger scale geographical structure, whereas subsequent spatial vectors represent smaller scales (Borcard & Legendre, 2002). Since MEMs describe how the sites are spatially distributed, the usage of such axes as variables also describes the spatial autocorrelation among plots. Our analysis resulted in 26 MEM axes, from which we selected subsets by fitting them against the response variables (i.e. alpha and beta diversities) via redundancy analyses (RDAs), retaining only those axes that were selected in a forward stepwise selection procedure (with function ordiR2step in package 'vegan ', Oksanen et al., 2019). From the selection procedure, up to 12 axes were selected for each of the response variables, with slight differences in the set of axes selected, but mostly representing geographical variation at larger scales (see maps of selected MEM axes in Appendix S2, Figure S2.3). To calculate the MEM, we used the function dbmem in R-package 'adespatial' (Dray et al., 2020).

| Traits and phylogenetic data
Trait data were obtained from the LEDA database (Kleyer et al., 2008).
We chose traits according to the leaf-height-seed scheme proposed by Westoby (1998) to capture essential ecological functions of plants, related to nutrient acquisition, competitive ability and establishment and dispersal respectively: specific leaf area (SLA), plant height (H) and seed mass (SM). SLA, defined as the ratio of total leaf area to total leaf dry mass, is positively correlated with the potential growth rate and represents the axis of resource exploitation ranging from more acquisitive to more conservative strategies (Wright et al., 2004). H is a whole-plant trait and it is related to dispersal ability, competition for light and expresses competitive ability, whereas SM is linked to both competitive ability during the seedling phase and dispersal (Díaz et al., 2016). When trait data were not available for a given species, we used the trait value of the closest taxonomic relative as an estimate, resulting in a complete species by trait matrix. Such a substitution was done for 48 of all 843 possible trait values (three traits across 281 species), that is, roughly 6%. Information of the phylogenetic relationships between species in the form of a phylogenetic tree was taken from Daphne (Durka & Michalski, 2012), a dated supertree of central and northwest European floras.

| Diversity measures
We calculated Rao's quadratic entropy index (Rao, 1982) for the FD and PD facets and Simpson's diversity for TD at alpha, beta and gamma levels, using the 'Rao' function provided by de Bello et al. (2010).
Within communities, Rao measures the diversity of species s in a sample based on the relative abundances p of the species, and the functional (i.e. trait) or phylogenetic distances d ij between the constituent For a given community, Rao is a generalization of Simpson diversity, as under the assumption of all species being functionally or phylogenetically unique (i.e. all d ij = 1) Rao be- As Simpson diversity, Rao can be calculated at different scales, corresponding to alpha (within communities), beta (among communities) and gamma (entire regional diversity) components, using an additive partitioning of gamma into alpha and beta (i.e. = − ), analogous to decomposing variance in an analysis of variance (ANOVA), as demonstrated by Pavoine and Dolédec (2005) and Ricotta (2005). An issue of many diversity indices is that their decomposition results in alpha and beta components that are not independent. Jost (2007) has suggested to remedy this issue via 'numbers equivalents', which for Simpson takes the form D = 1 ∕ (1 − D). The calculation of numbers equivalents is also applicable to Rao diversity, as demonstrated in detail by de Bello et al. (2010), and was adopted in our study. We used this approach (1) to decompose gamma diversity into alpha and beta components (in terms of proportional contribution), and (2) to calculate Simpson TD and Rao FD and PD values at the site level (alpha) and across the different sites (beta). For the site pairwise beta diversity, the same decomposition as described for the entire regional diversity above is used but applied to two plots at a time.
For FD we calculated Rao diversity indices for single functional traits (i.e. d ij representing differences between species for single traits), as well as for all three traits simultaneously (i.e. d ij representing average differences between species for all three traits), using the Gower distance for calculating distances.

| Statistical analyses
We analysed the variance of the three diversity facets and single trait diversities between sites explained by each set of predictors (climate, soil properties and geography) with RDA for alpha diversities, and distance-based redundancy analysis (dbRDA) for beta diversities. dbRDA allows RDAs when the response is a distance matrix, as in the case of pairwise beta diversities, through first performing a principal coordinate analysis on the distances (Legendre & Anderson, 1999). Climate, and soil and management were represented by the first two axes of their respective PCAs, whereas geography was represented by MEM axes (see above). Although the first two soil PCA axes jointly only accounted for 39%, including additional axes did not lead to better RDA model fits, which is why we did not include further axes. For each of the constrained RDA and dbRDA ordinations, the overall significance of the model, and the significance of the explanatory variables were assessed by Monte Carlo permutation tests, permuting the response variables 999 times.
Variance partitioning analyses were performed to examine the separate and combined effects of climate, geography and soil factors on alpha and beta components of the three diversity facets and single trait diversities. To evaluate the relationships between diversity facets with soil and climate, we tested the relationship between RDA results and individual environmental variables. For the alpha components, we tested the relationship between single diversity facets and the first two PCA axes of soil and climate variables, as well as geographical structure (via MEM selected axes), by fitting RDA models (function rda of R package 'vegan ', Oksanen et al., 2019). We used RDA instead of linear regression models to establish variable significance to maintain a uniform approach across alpha and beta level analyses. In any case, RDA with a univariate response as used here, is essentially analogous to multiple linear regression. However, the results were visualized by scatterplots and regression slopes from the linear regression model. For beta components, ordination diagrams indicating differences among sites according to environmental variables are given, in which also the individual MEM axes representing geographical space are included.

| Decomposition of diversity facets
Proportions of alpha and beta components differed among the diversity facets, with higher proportions of alpha in FD and PD (98.0% and 94.7%, respectively; Figure 2), while in TD contributions of alpha and beta diversity to overall gamma diversity were similar (46.0% and 54.0% for alpha and beta respectively).

| Variance partitioning of geography, climate and soil properties on diversity facets
Geography, environment (soil or climate) and the combination of geography and environment explained around 50% of the variation in both alpha and beta diversities, except for alpha FD, for which all the measured variables explained only 16% (Figure 3; Appendix S1, Table S1). Looking at single-trait FD, the environmental and spatial variables explained a greater amount of variation compared to multitraits FD (Appendix S1, Table S1). Spatially structured soil conditions (geography + soil) explained most of the variation in alpha TD (34%) and PD (34%), whereas alpha FD was mostly explained by geography (6%). Regarding single-trait FD, geography and spatially structured soil conditions were the variables that explained most of the variation for alpha diversity of all traits.
In addition to geography, spatially structured climate and soil conditions were also important depending on the facet considered.
Specifically, beta TD was secondarily explained by spatially structured climate (7%), whereas beta FD and PD were more explained by spatially structured soil conditions (12% and 28% respectively). Environment (either soil or climate) explained very little of both alpha and beta variation (Figure 3; Appendix S1, Table S1). Similarly, variation in beta single-trait diversity was mainly explained by geography and secondarily by spatially structured soil conditions (Appendix S1, Table S1).

| Alpha TD, FD and PD along soil and climatic gradients
Along soil gradients, alpha TD, FD and PD decreased from poor and basic soils towards more acidic and nutrient-rich sites ( Figure 4a; Appendix S1, Table S1). Moreover, alpha TD and PD increased when soil nitrogen and nitrate concentrations were higher (Figure 4b; Appendix S1, Table S1). None of the diversity facets were affected by climate at the alpha level (Appendix S1, Table S1).
Regarding single-trait diversity, H-FD and SM-FD decreased towards nutrient-rich and acidic soils, while SLA-FD increased. Along climate gradients we found much weaker but still significant increase in SLA-FD towards wetter communities characterized by higher temperatures, an increased H-FD in drier sites and SM-FD decrease towards lower temperatures ( Figure 5; Appendix S1, Table S1). None of the single-trait diversity indices were affected by nitrogen and nitrate concentrations (i.e. soil PCA2) at the alpha level (Appendix S1, Table S1).

F I G U R E 2
Proportion of alpha (in yellow) and beta (in blue) diversity from overall (gamma) diversity for the different diversity facets: Taxonomic diversity (TD), functional diversity (FD) and phylogenetic diversity (PD)

| Beta TD, FD and PD along space, soil and climatic gradients
Distance-based RDA diagrams showed how geography (represented by MEM axes), as well as climate and soil variables (PCA axes) affected species and lineages turnover among sites ( Figure 6, see statistics in Appendix S1, Table S1). Scattered and distant points in the diagram correspond to higher beta diversity between sites.
For both TD and PD facets, ordination of sites was mainly affected by large-scale geographical variables (MEM1-11). Conversely, depending on the facet considered, both spatially structured soil and climate variables were important. First, dbRDA analysis showed that both climate PCA axes (precipitation and temperature) were relatively important for beta TD, with more spread points in sites from western to north-western Europe (e.g. France and Belgium) related to temperature (climate PCA2) and dry conditions (climate PCA1). For beta FD, both climate and soil were relatively important, with higher beta FD related to temperature (climate PCA2-central Europe sites, e.g. German) and low soil pH and nutrient richness (soil PCA1-from western to north-western Europe, e.g. French and Belgium sites).
The ordination shows scattered points from central to northern Europe (e.g. Great Britain, France, Germany and Belgium), indicating higher beta FD for these sites, related to acidic and rich soils and both positive and negative values of climate PCA1 and PCA2. Finally, the F I G U R E 3 Proportion of explained variation in each diversity facet (FD, functional diversity; PD, phylogenetic diversity; TD, taxonomic diversity) for alpha and beta diversity explained by geography, climate, soil, combinations of two or all three (common), and residuals. In the bottom figures, the proportion of explained variation in single-trait diversity is shown (H-FD: Plant height diversity, SLA-FD: Specific leaf area diversity, SM-FD: Seed mass diversity) F I G U R E 4 Scatterplots indicating relationships between alpha diversity and the first two PCA axes of soil variables. (a) First PCA axis for soil variables, representing gradients of nutrients, pH, and management; (b) second PCA axis for soil variables, representing gradients in nitrogen and nitrate levels. Alpha diversity for the different facets (FD, functional diversity; PD, phylogenetic diversity; TD, taxonomic diversity) is scaled to the mean and centred around 0 to facilitate graph visualization. Note that relationships with single MEMs were omitted from the figure but reported in Table S1 ordination of sites in terms of beta PD was affected by both climatic and soil variables, with temperature being relatively more important (climate PCA2-central Europe, e.g. Germany). Results of single-trait diversity were included in Supporting Information (Appendix S1, Table S1; Appendix S2, Figure S2.4) since their contribution was marginal, showing similar patterns to the overall beta FD.

| DISCUSS ION
Although multiple faceted approaches at both alpha and beta levels are relatively common in ecology for different ecosystems and organisms (ants: Arnan et al., 2017;zooplankton: Gianuca et al., 2018;macroinvertebrates: Perez Rocha et al., 2018), such studies are surprisingly rare on European grasslands (but see Dainese et al., 2015).
To the best of our knowledge, our study is the first to integrate multiple facets at both alpha and beta levels to understand grasslands organization over a large geographical extent. Our results show that, by partitioning gamma diversity, Atlantic European grasslands of the

| Decomposition of overall diversity into beta and alpha components
Although species turnover is an important phenomenon among the studied communities, it is mediated by functionally similar and phylogenetically related species (i.e. functional and phylogenetic redundancy), since FD and PD turnover is low (Ricotta et al., 2020;Rosenfeld, 2002), and most of the overall variation in FD and PD can be found at the site (i.e. alpha) level. This overall low turnover indicates that climatic conditions are a strong regional environmental filter. Studying a particular European grassland type, not grassland vegetation in general, low FD turnover might be explained by similar climatic and edaphic conditions across our study sites, with species being filtered by similar traits to cope with these climatic conditions at local scales (see de Bello et al., 2009). On the other hand, low PD turnover might be related to the fact that grasslands are considered habitats of relatively young age, where only a few closely related lineages speciated during Paleogene and Neogene (Lososová et al., 2015). It is also worth noting that, although we are examining beta diversity within a specific grassland type, which F I G U R E 5 Scatterplots indicating relationships between functional alpha diversity and the first two PCA axes of soil and climate variables. (a) First PCA axis for soil variables, representing gradients in nutrients, pH, and management; (b) first PCA axis for climate variables, representing gradients in humidity; (c) second PCA axis for climate variables, representing gradients in temperature. Alpha diversity for the different traits (H-FD: Plant height diversity, SLA-FD: Specific leaf area diversity, SM-FD: Seed mass diversity) is scaled to the mean and centred around 0 to facilitate graph visualization. Note that relationships with single MEMs and soil PCA2 were omitted from the figure but reported in Table S1.3. Statistically significant single trait functional diversity trends (p < 0.05) are shown highlighted. MEM, Moran's Eigenvector Map; PCA, principal component analysis has been sampled based on 14 indicator species (Stevens, Duprè, Dorland, et al., 2011), it is unlikely that these species alone explain the overall low beta FD and PD, since the data comprise 281 species and span a large geographical range. Likewise, similar results of low beta FD despite high beta TD have been found in different studies and taxa, being mostly attributed to strong environmental filters selecting different species with similar strategies (Villéger et al., 2013;Wang et al., 2019). Moreover, the low contribution of beta to gamma does not indicate the absence of variation of beta diversity when we consider beta on a pairwise site by site basis. Indeed, the existing variation in beta diversity is explained by environmental variables, even to a degree that surmounts that with which variation in alpha diversity is explained (i.e. lower residual values of overall models of beta compared to alpha diversity, see Figure 3).

| Distinct importance of geographical, climatic and soil parameters
Soil parameters, despite being spatially structured at larger scales, are the most important filters at local scale (alpha diversity), confirming previous results found in the same grassland communities, suggesting that alpha TD is affected by large-scale gradients of soil properties (Stevens, Duprè, Gaudnik, et al., 2011). At larger scales (beta diversity), however, geography is the most important driver among all three facets, but spatially structured climate and soil conditions also contribute to beta diversity, depending on the facet considered.
Soil conditions have been recognized as important small-scale environmental filters in different ecosystems (Bernard-Verdier et al., 2013;Ceulemans et al., 2014;Helsen et al., 2014). Similarly, spatial distance and climatic conditions are expected to shape plant biodiversity at larger scales both because larger distances impose dispersal challenges to species and climate is an important environmental filter shaping regional species pools (Cornell & Harrison, 2014;de Bello et al., 2013). Interestingly, soil conditions affect only alpha TD and PD, not FD, reinforcing the importance of climatic filtering in European grasslands at larger scales. In other words, since beta FD is low and climate seems to filter similar traits at larger scales, consequently, soil conditions have a weaker role on alpha FD at smaller scales.
The contribution of geography driving beta diversity is often attributed to dispersal limitation, but it can also be explained by smaller scale stochastic processes such as ecological drift (Vellend et al., 2014), priority effects (Fukami, 2015) or even any unmeasured environmental variables (Peres-Neto & Legendre, 2010).
Therefore, stating the importance of dispersal processes based on spatial analysis needs to be done cautiously (Legendre et al., 2005;Peres-Neto & Legendre, 2010). However, as the examined dataset spans a large geographical area and only large-scale MEMs were selected, our results suggest that dispersal limitation and unmeasured environmental variables, not smaller scale stochastic processes, might be the most important drivers of beta diversity among all three facets. Indeed, different grassland studies have already shown the important role of dispersal limitation in driving grassland species diversity and composition (Leiva et al., 1997;Ozinga et al., 2005).
Geography, soil and climate explained around 50% of the overall variation in both alpha and beta diversity among all facets (except for alpha FD-see Appendix S1, Table S1). Such a result suggests that different unmeasured environmental variables, spatially structured or not, might be important to further explain the overall variation in alpha and beta TD, FD and PD across the study sites. For example, N deposition, and other anthropogenic effects, in general, are recognized as important drivers of vegetation change in grasslands Stevens, Thompson, et al., 2010) and could improve the power of detection in our analysis.

| Alpha diversity along gradients
The overall importance of soil at local scales can be attributed mostly to pH, soil nutrients and grazing, affecting all three facets.
The evidence of TD decreasing from basic towards acidic conditions is in line with studies based on analogous datasets, usually attributed to nitrogen deposition Stevens, Thompson, et al., 2010;Stevens et al., 2004). Moreover, the strong relationship between pH and soil nutrients, and the lower TD observed in more acidic and richer soils, confirms that soil acidification and eutrophication are driving biodiversity loss (Stevens, Thompson, et al., 2010). Nevertheless, we also expand previous knowledge of lower TD in richer and acidic soil in these grasslands Stevens, Duprè, Gaudnik, et al., 2011), indicating that both phylogenetic and FD (despite low effect size in the latter) are also lower in acidic soil. Therefore, high fertility, acidic and grazed sites in Atlantic European grasslands harbour fewer species, with lower diversity of traits and lineages. These results confirm patterns of decreasing FD with increasing nutrient availability found in the same study area (Helsen et al., 2014).
Lower diversity of all three facets in high fertility-grazing sites indicates that equalizing fitness, not limiting similarity or environmental filtering, is driving the community assembly of European grasslands at these smaller scales. Accordingly, although competition is often predicted to increase plant functional and/or PD, due to the exclusion of similar species that compete for the same resources (i.e. limiting similarity), we suggest that, towards fertile, acidic and grazed sites, stronger competitors (functionally and phylogenetically related) exclude weaker ones, decreasing the diversity of all three facets in those communities (Mayfield & Levine, 2010). Moreover, such decrease in all diversity facets might be connected to changes in proportion of growth forms composing the community, with an increasing proportion of grasses at the expense of forbs, as previously demonstrated in these grasslands Stevens, Duprè, et al., 2010).
Our results also show that alpha TD and PD are affected by the availability of different N forms in the soil. Accordingly, we found lower alpha diversity at soils richer in NH4+ and poorer in N and NO3−, confirming a higher availability of N in the form of ammonium in acidic soils (Pannek et al., 2015;Stevens, Thompson, et al., 2010) and the negative effects of NH4 on biodiversity (Pannek et al., 2015).
Studies have found either no effect (Stevens et al., 2004) or a negative effect of soil N on TD (Ceulemans et al., 2011) and FD (Helsen et al., 2014) in North-European Violion caninae habitats. Therefore, our results can confirm such findings, highlighting that soil N affects TD and FD facets but also PD, with different effects depending on the N form available in the soil.
Regarding the effect of land use on alpha diversity, since axis 1 of the soil PCA essentially ranges from grazed/mown to intensely grazed sites (Appendix S2, Figure S2.2), it must be noted how grazing had a negative effect on all facets (Figure 4). This pattern is in line with numerous studies highlighting pervasive effects of disturbance reducing grasslands TD, FD and PD (e.g. Dainese et al., 2015;de Bello et al., 2007;Egorov et al., 2014). As grazing by animals may greatly contribute to actual soil chemistry, and these two groups of variables are highly correlated, it is not possible for us to separate these effects. Likewise, an unbalanced proportion of mown (n = 35) and grazed (n = 118) sites is apparent across the study area, which in turn also reinforces the large proportion of overlap in explained variation from geography and soil portions.
By considering single trait variation along gradients, we show that most of the reduction in the overall FD found in high fertilitygrazing sites is due to height and seed mass, indicating that in those conditions, stronger competitors-likely taller and large-seeded species-outcompete weaker ones, as shown elsewhere (Dainese et al., 2012(Dainese et al., , 2015. Taller and large-seeded species may have competitive advantages in nutrient rich soils since they intercept light more efficiently and produce vigorous seedlings (e.g. Grime, 2006;Turnbull et al., 2004). Interestingly, the diversity in SLA increases in nutrient-rich communities likely due to different possible strategies that species adopt to exploit available resources, from conservative to acquisitive. These results reinforce the notion that single traits may respond differently during community assembly, being important to consider both multi-and single-trait variation (Bernard-Verdier et al., 2013;Spasojevic & Suding, 2012).
Although climate was not important for the multi-trait FD, single-trait variation suggests that towards wetter sites species become less diverse in terms of height and more diverse related to SLA. This result resembles the pattern found along fertilitygrazing gradients in our study, reinforcing the idea of equalizing fitness process in benign habitats (wetter in this case). We also show a decrease in SM-FD and SLA-FD towards lower temperatures, which indicates that climate acts as an important filter, selecting species with similar seed and leaf characteristics. Studies in mountain grasslands have also reported a decrease in seed mass and SLA-FD towards colder climates, suggesting that lower seed mass enables seeds to survive longer periods in the seed bank and low SLA is related to more conservative resource use strategies (Dainese et al., 2012;de Bello et al., 2013).  (Karger et al., 2015). Second, high beta diversity in less favourable climatic conditions (i.e. colder sites) might be related to the fact that only certain species are able to overcome either the abiotic or biotic filtering operating in those conditions (de Bello et al., 2013). Therefore, the importance of climate driving taxonomic, functional and phylogenetic composition at larger scales in European grasslands could be attributed to both competition and environmental filtering.

| Beta diversity along gradients
Relatively high beta FD and PD in poor/basic soil is correlated with higher alpha FD and PD in the same soil conditions. This result could be explained in a probabilistic way, indicating that, as the overall beta FD and PD is low, the chances of finding variation in traits and lineages are higher in sites where the diversity of traits and lineages is higher. Interestingly, on the opposite side of the soil gradient (acidic/rich sites) FD and PD turnover is low and equalizing fitness is an important process driving the community assembly at smaller scales. This result suggests that the equalizing fitness processes operating in European grasslands are being mediated by strong competitors with a similar set of traits and evolutionary history. The combination of low overall beta FD and PD and similar strong competitors may indicate that accelerating eutrophication in European grasslands, due to high levels of N deposition Stevens, Thompson, et al., 2010), might increase biodiversity loss and homogenization among all facets not just at local scales but also at large scales.

| CON CLUS IONS
Our results highlight the importance of integrating multiple biodiversity facets at both alpha and beta levels, as well as considering both multi-and single-trait variation along environmental gradients, to advance theoretical and practical issues in plant community assembly studies. Overall, the community assembly on Atlantic European grasslands is mediated by functionally and phylogenetically similar species, driven at larger scales by stochastic processes (i.e. dispersal limitation) and climatic filters, and at smaller scales by soil filters (pH and nutrients) and equalizing fitness processes at nutrient-rich sites. Moreover, from a conservation perspective, our results suggest that accelerated acidification and eutrophication of European grasslands impose threats not only to TD, as previously reported Stevens, Duprè, et al., 2010;Stevens, Thompson, et al., 2010), but also to both alpha and beta FD and PD.
Therefore, to maintain long-term grasslands stability and functioning, biodiversity conservation of Violion caninae communities, and grasslands in general, should be considered at multiple biodiversity facets and broader scales.

ACK N OWLED G EM ENTS
We thank the reviewers for their input in improving this manuscript and Lorella Dell'Olmo (University of Florence) for the help in the production of maps. No permits were required for this study.

CO N FLI C T O F I NTE R E S T
The authors have no conflict of interest to declare.

DATA AVA I L A B I L I T Y S TAT E M E N T
The vegetation, soil, and management data were retrieved from Stevens, Duprè, Dorland, et al. (2011) Author contributions: All the authors conceived the ideas; MM, DPFT, MT and LG analysed the data; MM, DPFT, MT, KK and LG led the writing and JH revised and critically improved the manuscript.

AUTH O R B I O G R A PH I E S
Michele Mugnai is a PhD candidate; he is interested in community ecology, with special focus on plant functional diversity.
Diego P. F. Trindade is a PhD candidate interested in understanding how plant communities are organized across space and time.
Mélanie Thierry is a researcher in community ecology.
Krishan Kaushik is a PhD candidate, he works in temperate plant communities involving invasion and functional trait perspectives.

Jan
Hrček is a lab leader with focus on experimental community ecology & evolution.
Lars Götzenberger is a lecturer and researcher in trait-based plant ecology.

S U PP O RTI N G I N FO R M ATI O N
Additional supporting information may be found in the online version of the article at the publisher's website.