Positive Effects of Legumes on Soil Organic Carbon Stocks Disappear at High Legume Proportions Across Natural Grasslands in the Pyrenees

Soil is the largest terrestrial carbon pool, making it crucial for climate change mitigation. Soil organic carbon (SOC) is suggested to depend on biodiversity components, but much evidence comes from diversity-function experiments. To disentangle the relationships of plant guild diversity with SOC storage (kg m−2) at broad spatial scales, we applied diversity-interaction models to a regional grassland database (n = 96) including wide environmental conditions and management regimes. The questions were: (1) Are the effects of plant guilds on SOC stocks in natural grasslands consistent with those found in experimental systems? (2) Are plant guild effects on SOC stocks independent of each other or do they show interactive—synergistic or antagonistic—effects? (3) Do environmental variables, including abiotic and management, modify guild effects on SOC stocks? Among our most novel results we found, legume effects on grassland SOC vary depending on legume proportion consistently across broad spatial scales. SOC increased with legume proportion up to 7–17%, then decreased. Additionally, these effects were strengthened when grasses and forbs were codominant. Grazing intensity modulated grass proportion effects on SOC, being maximum at relatively high intensities. Interpreting our results in terms of existing contrasted ecological theories, we confirmed at broad spatial scales and under wide-ranging environmental conditions the positive effects of plant guild diversity on SOC, and we showed how legumes exert a keystone effect on SOC in natural grasslands, probably related to their ability to fix inorganic N. Niche complementarity effects were illustrated when codominance of forbs and grasses at optimum legume proportions boosted SOC storage, whereas grass dominance increased SOC stocks at medium–high grazing intensities. These findings can facilitate the preparation of regional and local strategies to ameliorate the soil capacity to absorb carbon.

antagonistic-effects? (3) Do environmental variables, including abiotic and management, modify guild effects on SOC stocks? Among our most novel results we found, legume effects on grassland SOC vary depending on legume proportion consistently across broad spatial scales. SOC increased with legume proportion up to 7-17%, then decreased. Additionally, these effects were strengthened when grasses and forbs were codominant. Grazing intensity modulated grass proportion effects on SOC, being maximum at relatively high intensities. Interpreting our results in terms of existing contrasted ecological theories, we confirmed at broad spatial scales and under wide-ranging environmental conditions the positive effects of plant guild diversity on SOC, and we showed how legumes exert a keystone effect on SOC in natural grasslands, probably related to their ability to fix inorganic N. Niche complementarity effects were illustrated when codominance of forbs and grasses at optimum legume proportions boosted SOC storage, whereas grass dominance increased SOC stocks at medium-high grazing intensities. These findings can facilitate the preparation of regional and local strategies to ameliorate the soil capacity to absorb carbon.

INTRODUCTION
Soil carbon is highly relevant for climate change mitigation (Canadell and others 2007) because it constitutes the largest carbon pool in terrestrial ecosystems (Batjes 1996). According to current conceptual models, soil carbon sequestration is expected to depend strongly on biodiversity components (Hooper and others 2005;Díaz and others 2007;Hector and Bagchi 2007;Kirwan and others 2009;Suter and others 2015;Connolly and others 2018). However, most evidence comes from diversity-function experiments, which might suggest important underlying mechanisms (Fornara and Tilman 2008; Prommer and others 2020), but do not encompass the whole complexity encountered in real ecosystems. In particular, experimental studies tend to neglect interactions and to reflect only short-term processes; thus, extrapolating their results for real systems could lead to incomplete or misleading conclusions (Yuan and others 2017;Rillig and others 2019). In contrast, evidence about a positive association between soil carbon storage and plant diversity in observational, broad-scale studies is scarce despite some recent attempts (Hollingsworth and others 2008;Manning and others 2015;Wardle 2016;Chen and others 2018).
Although the role of every plant species in the ecosystem may be unique (Wayne and Bazzaz 1991), grouping them into plant guilds provides good summary representations for biodiversityfunction analysis (Sebastià 2007). Plant functional types (PFT) are suites of species that share ecological characteristics and play similar roles in ecosystems (Steneck, 2001;Blondel, 2003). When defined in terms of resource use, they are often called guilds (see Sebastià 2007;and Root 1967 for animals).
Among other biodiversity approaches, plant guilds have been found to be highly relevant for explaining soil carbon content in experiments (Fornara and Tilman 2008;Lange and others 2015). Considering grassland ecosystems, grasses are efficient in terms of light capture because their leaves are at vertical angles (Sebastià 2007;Zhou and others 2016), while the architecture of their roots makes them efficient capturing soil N, among other nutrients. Legumes can have access to symbiotically-fixed atmospheric N and have recently been found to be more active in terms of CO 2 exchange per biomass unit compared with other guilds (Ibanez and others 2020). Nonlegume forbs present a variety of ports and leaf angles, can store high amounts of nutrients in their roots and cannot fix atmospheric N (Schmid and others 2004;Lipowsky and others 2015;Herz and others 2017). Finally, sedges share multiple morphological features with grasses, but often present some differences in their photosynthetic tissues and higher root to shoot ratios (Mou and others 2018). The consideration of sedges as a separate group from grasses was recommended by Onipchenko and others (2012), because of their high nutrient uptake capacity over grasses in that study.
Plant guild effects on soil carbon storage need to be studied at broad scale to understand how they work under different broad-scale abiotic conditions, as is being done for plant taxonomical and functional diversity (Manning and others 2015;Carol Adair and others 2018). Ecologists and modelers need this information to validate conceptual paradigms and generate new hypotheses contributing to the refinement of global mechanistic models. Land managers and policymakers need it to establish priorities for conservation objectives. Grasslands are ideal to study soil organic carbon (SOC)-biodiversity interactions, due to their wide distribution, high SOC stability and biodiversity compared to other ecosystems (Suttie and others 2005;Poeplau and Don 2013;Dengler and others 2014). Taking into consideration that previous studies, carried out at narrower scales, found a crucial but contrasting role of legumes among the plant guilds (e.g., Fornara and Tilman 2008;Lange and others 2014), we aim to answer the following questions: 1. Are the effects of plant guilds on SOC stocks in natural grasslands consistent with those found in experimental systems? 2. Are plant guild effects on SOC stocks independent of each other or do they show interactive-synergistic or antagonistic-effects? To answer these questions, we built statistic models of SOC stocks distribution using an extensive database generated from a survey of 96 natural grasslands in the Pyrenees. The database includes a variety of climates; different landscape positions; and a range of grazing management regimes. This environmental heterogeneity is also reflected in the wide variability of plant guild composition (Sebastià 2007;Gó mez 2008). We applied the diversity-interaction modeling (DIM) approach (Kirwan and others 2009), which allowed us to separate the identity and interaction effects of plant guilds. This approach consists of testing a series of hypotheses, from the absence of biodiversity effects on SOC to the existence of interaction effects involving both plant guilds and environmental drivers (Kirwan and others 2009;Connolly and others 2013).

Study Sites
The set of data used in this study was extracted from the PASTUS Database (http://ecofun.ctfc.cat/ ?p=3538), compiled by the Laboratory of Functional Ecology and Global Change (ECOFUN) of the Forest Sciences Centre of Catalonia (CTFC) and the University of Lleida (UdL). We sourced a wealth of data of 96 grassland patches distributed across the Pyrenees ( Figure S1), including topographical, climate, soil, herbage and management variables. The elaboration of the PASTUS Database concerning this study is summarized in Figure S2. The sampled area encompasses a wide variety of temperate and cold-temperate climates with different precipitation conditions, depending on altitude and geographical location, from Mediterranean to continental and Boreo-Alpine environments (de Lamo and Sebastià 2006;others 2018, 2020). Almost all of the plant species in the grasslands from the PASTUS Database are perennial (Sebastiá and Sebastià 2004), and both plants and environments are highly heterogeneous (Rodríguez and others 2018).

Sampling Design
Sampling in the PASTUS Database was designed according to a stratified random scheme, where samples were selected at random within strata. This process was undertaken using the ArcMap 10 software (ESRI, Redlands, CA, USA). The basis for randomization was the map of habitats of Catalonia 1:50 000 (Carreras and Diego 2006) for the eastern and central sectors of the Pyrenees; the map of habitats of Madres-Coronat 1:10,000 (Penin 1997) for the northeastern sector; and the land use map of Navarre 1:25,000 (Gobierno de Navarra 2003) for the western sector. Four variables were initially considered for sampling stratification within each sector: altitude (< 1800; 1800-2300; > 2300 m), slope (0-20; 20-30; > 30°), macrotopography (mountain top/south-facing slope; valley bottom/ north-facing slope) and grazer type (sheep; cattle; mixed). Accordingly, we determined a set of homogeneous grassland patches by crossing the stratification variable layers. Grassland patches were then listed by type and arranged within each list randomly to determine sampling priority. At least one to two replicates of each patch type, defined by the stratification variables, were sampled.
In each sampled grassland patch, a 10 9 10 m 2 plot was systematically placed in the middle of each homogeneous grassland patch, including a particular plant community. We collected soil and vegetation samples and assessed environmental variables within each 100 m 2 plot (see others 2018, 2020), for additional details about sampling design). Local environmental variables were assessed inside the 100 m 2 plots.

Climate, Topographical and Grazing Drivers
All the variables included in this study are listed in Table S1. Climate variables were determined from Worldclim 2.0 (Fick and Hijmans 2017). Based on previous modeling of SOC in grasslands in the Pyrenees (Rodríguez and others 2020), we selected mean annual precipitation (MAP) and the intraannual difference of temperature (Temperature Seasonality Index of Sebastià (TSIS); (Debouk and others 2020; Rodríguez and others 2020)), calculated as the difference between mean summer temperature and mean annual temperature.
Topographical variables included slope and macrotopography. Slope was determined in the field by clinometer. Macrotopography included exposed (south-and east-facing slopes, mountain tops) and protected (north-and west-facing slopes, valley bottoms) positions.
Regarding grazing management variables, detailed surveys were carried out among farmers, shepherds and land managers, to complete the initial general information used in the sampling design. We used these data to determine grazing intensity by estimating livestock stocking rates measured as livestock units/ha (LU/ha). We treated it as a semiquantitative variable with three categories (Sebastià and others 2008): low (1: < 0.2 LU/ha), moderate (2: 0.2-0.4 LU/ha) and highmoderate (3: up to 0.4 LU/ha); a few samples corresponded to abandoned grasslands (0 LU/ha).

Soil Environmental Factors and SOC Stocks
In each plot, one soil core was extracted with a 5 9 5 cm probe at 0-20 cm soil depth. To determine bulk density, we air-dried and weighed the soil samples: We then sieved each sample to 2 mm to separate stones and gravels from the fine earth fraction. The fine fraction was sent to the laboratory to determine soil texture, pH and soil organic carbon (SOC) stocks. Textural classes in the top 10cm soil layer were determined by the Bouyoucos method (Bouyoucos 1936), and soil pH in the top 10-cm soil layer was measured in water (Solly and others 2014). Total carbon (C) concentration of the fine earth was determined by the elemental autoanalyzer. Soil organic carbon stocks (SOC) in the upper 20-cm soil layer were then estimated taking into account the organic C concentration in the sample and its bulk density, and subtracting the coarse particle (> 2 mm) content, following García-Pausas and others (2007). Hence, we obtained SOC stocks in kg m 2 at a fixed depth of 20 cm. See Rodríguez and others (2020) for further details about SOC sampling and determination.

Plant Guild Proportions
Plots of 10 9 10 m 2 were established in the middle of homogeneous grassland patches holding a given plant community (Canals and Sebastià 2000;Sebastiá 2004). Aboveground biomass was estimated from herbage harvested in four 50 9 50 cm quadrats placed in a 2 9 2 m subplot within the 100 m 2 plot. Plant guild biomass was determined in one of the four quadrats per plot by hand separation. Four guilds were sorted: C3 grasses, sedges, legumes and non-legume forbs (the latter including some subshrubs), following Sebastià (2007). Guild biomass was then oven-dried at 60ºC to constant weight.

Diversity-Interaction Models
To disentangle plant guild effects on SOC, we built a series of models following the diversity-interaction model (DIM) approach (Kirwan and others 2009). We log-transformed SOC (logSOC) using the natural logarithm to prevent a breach of the normality assumption by the residuals of the models ( Figure S3). DIM models were fitted in a four-step process based on equations from Kirwan and others (2009) and Connolly and others (2013). First, we defined a null model; Model 0: where logSOC (y) depends on a series of terms including environmental abiotic variables and interactions between them (X), each one affecting y according to its coefficient b. The remaining variance of y is provided by model error . According to this null model, there are no effects of plant diversity, including plant guild proportions, on SOC. We decided to use the geophysical model presented in Rodríguez and others (2020), using data from the same database, as this null model. Although the set of data employed for fitting the two models was not exactly the same because plant guild variables were not measured in all the PAS-TUS Database samples, this model properly explained a significant amount of the total SOC variance (R 2 adj = 0.37) when including only the sites selected for the present study. The model terms of this model (Table 1) were already explained in Rodríguez and others (2020), allowing us to focus on plant guild effects in this work. Additionally, models obtained with the same procedure using only the sites in this work (Appendix) included only MAP, TSIS and clay content, whereas the model from Rodriguez and others (Rodríguez and others 2020) finally used here also included topography and grazing management variables. This is an advantage since the variables of the null model will be used to test interactions with plant guilds in Model 4.
Secondly, we fitted Model 1: which simply consists of the addition of each plant guild proportion P and the estimation of its effect on y (b) to the null model. This model assumes that plant guilds affect SOC exclusively by their single (identity) effects, without interactions. Note that these models are similar to a conventional linear model (Legendre and Legendre 1998) but removing the intercept to allow the inclusion of all the plant guild proportions (Connolly and others 2013).
Model 2 (evenness model) assumes that there are interaction effects between plant guilds and that they are the same for all the pairwise plant guild combinations: is a measure of the relative abundance of plant guild proportions and is maximum when guilds are all equally represented others 2009 and2007). d AV is the single common interaction coefficient. In contrast, Model 3 (separate pairwise interactions) hypothesizes that effects of pairwise interactions between plant guilds may be different from each other and dependent on the relative abundance of each plant guild: Thus, the term d ij indicates the strength of the interaction between the plant guilds i and j, which is modulated by the relative abundance of both groups (P i and P j ).
Model 4 adds the possibility of including pairwise interactions between plant guild proportions and abiotic environmental variables: We selected P i X candidate terms with a semiautomatic approach, following Rodríguez and others (2020). We used the glmulti function from the glmulti package (Calcagno 2015) with the genetic algorithm option. We included the explanatory variables of Model 0 (TSIS, MAP, macrotopography, slope, clay content and grazing intensity) and the three most abundant plant guilds (grass, forb and legume) as explanatory variables, and the residuals of the geophysical model from Rodríguez and others (2020) as a response variable (i.e., Model 0 with intercept, since glmulti function requires it). This semiautomatic process began by obtaining a best subset of models using the corrected Akaike information criterion (AICc), appropriate when n/k is less than 40, n being the sample size and k the number of parameters in the most complex model (Symonds and Moussalli 2011). We selected the best model and its equivalents obtained by the genetic algorithm, which were those within two Akaike information criterion-corrected (DAICc) values of the best model, as a DAICc < 2 indicates that the candidate model is almost as good as the best model (Burnham and Anderson 2002). MAP: mean annual precipitation; TSIS: temperature seasonality; slope: terrain slope; exposed: exposed position according to macrotopography; clay: clay content; low and moderate intensity: low and moderate grazing intensity.
For this set of models, we built averaged models using the MUMIn package (Barton 2015). We calculated partial standardized coefficients, obtained by multiplying the unstandardized coefficient in the model by the partial standard deviation of the variable, which is a function of the standard deviation of the variable, the sample size, the number of variables in the model and the variance inflation factor of the variable (Barton 2015). We selected all the variables with significant effects (alone or in interaction with each other) in the conditional average model, which was preferred over the full average model because we wanted to avoid excessive shrinkage effects at this moment of the modeling procedure (Grueber and others 2011). To finally obtain Model 4, we added these terms to Model 3, and removed those which were not significant according to an F-test by the anova function in R, comparing the model with one lacking the corresponding term (de Vries and others 2012).
Model 5 inserts the hypothesis of a triple interaction effect among the plant guilds: l ijk being a coefficient modulated by three of the four plant guilds (P i P j P k ). Finally, Model 6 includes a term (h) that represents the shape of the interaction relationship, being a straight line when h = 1. This model was fitted with the nls R function, using the coefficients of Model 5 as starting values (Connolly and others 2013).
Once the models 0-6 were built, we compared them sequentially by F-test using the ''anova'' R function and the likelihood ratio test (''lrtest'') function in the lmtest R package (Hothorn and others 2019). These analyses facilitated the test of whether the terms added in a model provided significant information with respect to the previous model. Additionally, we calculated the AIC difference of models 1-5 with Model 0 (DAIC) to ac-count for the improvement provided for each model, including plant guilds compared to the null model. Among models 0-6, we selected the best model considering two criteria. First, as models 0-6 are sequentially and hierarchically more complex, the best model should significantly improve its precedent model. Hence, its corresponding F and likelihood ratio tests were required to be significant. Second, the best model could not be considered equivalent to Model 0. To account for this statement, DAIC was required to be higher than 2 (Burnham and Anderson 2002).
We used the emmeans package (Lenth and others 2019) to estimate SOC stocks under certain conditions according to the selected model. Firstly, in order to show a complete perspective of the effects of plant guilds, we calculated the predicted effects of the model for all possible combinations of grasses, forbs and legumes, fixing the rest of the variables at their mean value. Sedges were fixed at 0% because of their low effects on SOC, since they were not involved in any significant interaction (Table 2) and given their low mean proportion (2%). We built a ternary plot with the ggtern package (Hamilton and Ferry 2018) to show predicted SOC variation across plant guild proportions. Secondly, to discuss the effects of legumes at different grass and forb proportions, we used the same methodology to plot SOC variation across a legume proportion gradient under: (1) codominance scenario where grasses and forbs accounted for equal proportions (50:50 of the non-legume proportion); (2) a forb-dominance scenario where forb/grass proportions were 80:20 of the non-legume proportion; and (3) a grass-dominance scenario where forb/grass proportion was 20:80 of the non-legume proportion. Thirdly, to study the effect of grasses at different grazing intensities, we drew SOC variation across a grass proportion gradient at the three grazing intensity levels under a forb:legume codominance scenario. We complemented the information provided with the last two plots with the changes in the slope produced along the corresponding legume and grass proportion gradients. We calculated slope at each point as Slope ¼ y 2 À y 1 x 2 À x 1 y 2 and x 2 being the values corresponding to each integer value of x, and y 1 and y 1 the values of the point x 2 -1 (Nagle and others 2017).

RESULTS
The best model explaining SOC stocks included the triple interaction between grasses, legumes and forbs, in addition to plant guild proportion and pairwise interactions between grasses, legumes and forbs. Furthermore, the interaction between grass proportion and grazing intensity also entered the best model, with a DAIC of -15.56 (Model 5; Tables 2 and 3). Testing the possibility of nonlinearity of interaction effects of Model 5 (that is, Model 6) did not improve the model (Chi-square = 12.7; p value = 0.44; Table 2). The effects of sedges were independent of other plant guilds and environmental factors, yielding extremely low SOC stocks (Table 3). Considering these results and the scarcity of sedges in our dataset, we focus on the other three guilds. SOC stocks were maximized at balanced grass and forb mixtures (40% to 60% grasses and 20% to 45% forbs, approximately) combined with moderate legume proportion (10% to 30%, approximately), according to Model 5 (Table 3 and Figure 1). Legumes enhanced SOC storage when mixed with either grasses or forbs, but only at low to moderate legume proportions (Figure 2). SOC storage increased with legume proportion, up to 7-17% legumes, depending on neighbors (Figure 2). At legume proportions above 7-17%, SOCS decreased with increased legume proportion (Figure 2). This enhancement effect was maximum when grasses and forbs were at an equivalent proportion, whereas it was moderated in the grassdominance scenario and strongly reduced in the forb-dominance scenario.
On the other hand, the positive effects of grass proportion on SOC stocks were strongly conditioned by grazing intensity, modifying the patterns of the SOC-plant guild relationship (Figure 3; Table 3). At moderate grazing intensities ( Figure 3B), SOC stocks reached their maximum values at 7-17% legume proportions and 40-60% of grass and forb proportions. In contrast, at high-moderate intensities (Figure 3), SOC increased with grass proportions, reaching maximum values in grasslands with 100% of grasses. Moreover, the SOC grasses relationship was always positive at this grazing regime ( Figure S4). Finally, the SOC pattern at low grazing intensities was intermediate between the previous ones (Figures 3 and S4

Consistency Between Plant Guild Effects on SOC Stocks in Natural Grasslands and Experimental Systems
SOC storage represents a complex equilibrium between primary production (inputs) and organic matter decomposition (outputs) that depends on multiple environmental factors, including climate, soil texture and nutrients, and land management (Jenny 1941;Schlesinger 1977;Jackson and others 2017). Our results demonstrate that some guilds at low proportions can have positive effects on SOC, but also that those enhancement effects can disappear or shift when the proportions among guilds are modified. In our DIM (Table 1 and Figure 1), the critical role of legumes on SOC is not surprising, since their symbiosis with rhizobia bacteria allows them to fix N 2 from the atmosphere (McGrath and others 2014), having large effects on N availability and supply (Zanetti and others 1997;Spehn and  -7 9 10 -5 2.4 9 10 -4 -0.28 0.78 Grasses x Legumes 1.1 9 10 -4 3.4 9 10 -4 -0.28 0.75 Grasses x Forbs -1.3 9 10 -4 7.6 9 10 -5 -1.77 0.08 Sedges x Legumes 5 9 10 -5 6.2 9 10 -4 0.9 0.932 Sedges x Forbs -1.1 9 10 -4 2.7 9 10 -4 -0.39 0.696 Legumes x Forbs -8 9 10 -5 4.8 9 10 -4 -0.17 0.868 Grasses x Forbs x Legumes 2.3 9 10 -5 9 9 10 -6 2.42 0. Grasses: proportion of grasses; Sedges: proportion of Sedges; Legumes: proportion of legumes; Forbs: proportion of forbs; MAP: mean annual precipitation; TSIS: temperature seasonality; slope: terrain slope; exposed: exposed position according to macrotopography; clay: clay content; low and moderate intensity: low and moderate grazing intensity. These effects on soil N may alter SOC inputs and outputs among other ecosystem processes (Hector and others 1999;Fornara and Tilman 2008). However, what is original in our work is the changing role of legumes on SOC storage depending on plant diversity and plant guild interactions. Former studies have reported positive, neutral and negative effects of legumes on SOC (Steinbeiss and others 2008;Lange and others 2015;Wu and others 2017). Results from our DIM are novel because they suggest that those effects could be dependent on legume biomass proportion, as found in the case of Pyrenean natural grasslands. others (2009, 2011) found similar responses of plant N and yield to legume proportions in an experimental work, which may be behind the SOC-legume relationship we found. At low proportions (0-17%), an increase in legumes noticeably enhanced the positive effect on SOC (Figures 1 and 2), which could be classified as a keystone effect (Mills and Doak 1993). Positive effects of legumes on SOC may have three different non-excluding explanations according to Zhao and others (2014). Firstly, legumes would promote an increase in primary productivity of plant communities through increased N availability, which would lead to an in-crease in SOC (Wu and others 2017). Secondly, positive effects of legumes on SOC can be attributed to low C/N ratios of legume residues (i.e., litter, root exudates…), more similar to soil microorganisms and soil organic matter than other plant residues (Jensen and others 2012). Substrates with low C/N ratios can reduce microbial N acquisition and increase their carbon use efficiency, facilitating humification processes and plant residue decomposition into soil organic matter (Spohn and others 2016). Thirdly, N inputs through legumes may inhibit the production of oxidative enzymes to degrade the more recalcitrant compounds, leading to reduced ecosystem CO 2 emissions and C losses (De Deyn and others 2011; Spohn and others 2016).
On the other hand, in our study, when legume proportions were high, they had a negative effect on SOC according to our model (Figures 1 and 2). Negative effects of legumes on SOC stocks have been attributed to a decrease in community root biomass (Lange and others 2015). Recent work by Prommer and others (2020) suggests that the allocation of carbon to roots is less necessary at the community level when legumes are present, due to the increase in available N they fix. Furthermore, Henneron and others (2019) found that the resource acquisitive strategy of legumes, with high photosynthetic activity and capacity and high root metabolic activity, exudation and death, may enhance soil microbial activity, depleting SOC stocks (i.e., rhizosphere priming effect; Kuzyakov, 2002). This explanation is also suggested by the findings in Ibañ ez and others (2020), where high efficiency of legumes in C capture does not translate into high legume biomass in the grassland. We also suggest that, at some point, pathways for positive effects of legumes on SOC may be reversed. For instance, due to legumes strategy of producing a more nutrient-rich and short-lived biomass than other plant guilds (Craine and others 2002), high legume proportions may involve less total biomass in the long term, and very low litter C/N ratios could lead to higher mineralization rates and less carbon storage (Orwin and others 2010). Despite the fact that labile plant litter components can be a source of microbial-transformed products, and hence, of soil organic matter (Jensen and others 2012; Cotrufo and others 2013), some studies suggest that high-quality litter may have lower incorporation rates into the most stable soil organic matter than low-quality litter when the plant carbon inputs are high enough (Castellano and others 2015; Sollenberger and others 2019). Indeed, this fits in with findings by Ibañ ez and others (2020) in grasslands in the Pyrenees, where higher legume proportions are associated with lower yield in spite of high legume photosynthetic efficiency. Moreover, at high legume proportions, the additional N supply provided by N 2 fixation may be inhibited because of the reduction in plant competition for soil N (Nyfeler and others 2011); hence, the mechanisms which may enhance SOC stocks at low legume proportions could be inhibited at high legume proportions.

Interactions Between Plant Guilds Shaping SOC Stocks
Importantly, this enhancement effect on legumes over SOC stocks depended on the relative proportions of grasses and forbs, being maximum in the codominance scenario and minimum in the forb dominance scenario, with the grass-dominance scenario in between ( Figure 2). Maximum enhancement effects in the codominance scenario suggest a positive diversity effect on SOC stocks. According to experimental evidence, diversity produces spatial (different root depths) and temporal (different growing seasons) niche partitioning (Van Ruijven and Berendse 2005; Mueller and others 2013), which would allow the maximization of resource use (Blair and others 2014), leading to an increased plant productivity (Mason and others 2020) and microbial activity (Lange and others 2015) that would enhance SOC stocks (Chen and others 2018). Importantly, if legume proportions are close to 0, the variation in SOC stocks due to changes in grass and forb proportions is highly diminished (Figure 1). Hence, besides the keystone effect of legumes, our model suggests a niche complementarity effect (Diaz and Cabido 2001) represented by the triple interaction term (Table 3). The fact that grasses have slightly better performance on SOC stocks than non-legume forbs could be related to positive synergistic effects between grasses and legumes concerning N fixation and yield (Kirwan and others 2007;Nyfeler and others 2011;Rasmussen and others 2012;Schipanski and Drinkwater 2012;Ribas and others 2015;Suter and others 2015). Pirhofer-Walzl and others (2012) suggested that grasses could be especially good receivers of legume-derived N due to their fibrous root systems, which provide grasses with larger root surface and superficiality. Additionally, N 2 fixation activity of legumes could be promoted by an increased demand for soil N on the whole ecosystem, as a consequence of an enhancement of grass root systems driven by the supply of atmospheric N 2 fixed by legumes and transferred by litter, dead roots and exudates (Nyfeler and others 2011). That enhancement of grass root systems would also allow grasses to acquire more N from non-symbiotic sources (Nyfeler and others 2011;Suter and others 2015).

Interaction Effects of Grasses and Management Intensity on SOC Stocks
In addition to the interaction effects with legumes and forbs, grasses were the only plant guild which interacted with an environmental factor to shape SOC stocks. Interactions between vegetation and other variables affecting SOC stocks have been anticipated (Chang and others 2018;Chen and others 2018;Yuan and Jiang 2021). However, a recent study in the Northern Iberian Peninsula reported that only soil microbial nitrogen was affected by interactions between plant guilds and environmental factors, among a series of soil activity and fertility parameters involved in the C and N cycles (Debouk and others 2020). The effects of grass proportions on SOC stocks varied depending on grazing intensity level (Figure 2), a factor known to affect SOC through different paths, such as regulating the quantity and quality of organic matter returned to soil (Bardgett and Wardle 2003) or affecting soil respiration and nutrients by animal trampling or soil microbiota alteration (Lu and others 2017). However, positive effects of grazing on SOC stocks are often related to the promotion of aboveground and root biomass production (Franzluebbers and others 2000; Zeng and others 2015), and therefore this interaction effect could be related to the higher regrowth rates and productivity of grasses under moderate grazing pressures, when compared to the other groups Ganjurjav and others 2019). Our model suggests that only under relatively high grazing intensities (in our natural grasslands all grazing intensities are actually low; see methods) can grass dominance (< 60%) overcome the SOC decrease caused by the loss of the equilibrium between grass, forb and legume proportions (Figure 1).

Implication of Results in the Current Ecological Framework
Most of the previous studies addressing plant guild effects on SOC of grasslands were carried out at local scales and/or employing experimental assemblages (Fornara and Tilman 2008;Prommer and others 2020), not natural ecosystems as considered here. Moreover, most of these studies only considered the effect of plant guild richness and the presence or absence of the different plant guilds on SOC (Lange and others 2015; Wu and others 2017), although guild proportion effects have been described for other ecosystem functions, including yield (Kirwan and others 2007;Finn and others 2013;Ribas and others 2015). In contrast, what our model suggests is that: plant guild effects vary depending on their mass proportion in plant communities; some guilds, even in small proportions, can greatly enhance ecosystem function; those effects can be reversed at increasing proportions; and those observed effects stand over a wide range of grasslands and environmental conditions at regional scale. Hence, in addition to direct effects, shifts in plant guilds will have clear consequences on the rest of the biota (i.e., other plants, microbes), triggering a set of cascading effects in the ecosystem (Loranger-Merciris and others 2006; Cornwell and others 2008;De Deyn and others 2008). We postulate that legume proportion determines the predominance of some processes over others, leading to SOC accumulation at moderate legume proportions and to SOC depletion at high legume proportions.
Our results are also relevant for functional ecology, since they illustrate the power and usefulness of plant guilds to study keystone effects and interactions in ecosystems. In the last few decades, PFTs or plant guild approaches have been described as inferior methods in comparison with continuous trait indexes (Mason and others 2005;Funk and others 2017). However, these last methods focus on testing concrete hypotheses, like the mass ratio (Grime 1998;Díaz and others 2007) or niche complementarity (Villé ger and others 2008; de Bello and others 2016) hypotheses. In the present study, we detected other critical functional effects of plant guilds on SOC storage, like keystone legume effects (Spehn and others 2002) or guild interaction effects (Fry and others 2014). Additionally, our work demonstrated that plant guilds are also valid to detect niche complementarity or diversity effects (Diaz and Cabido 2001), because intermediate legume, grass and forb proportions enhanced SOC stocks, as well as mass ratio or dominance effects (Grime 1998;Díaz and others 2007), since SOC stocks reached maximum values under grass dominance at relatively high grazing intensities. Thus, our results suggest that particular environmental conditions could modify guild (biodiversity) interaction effects on ecosystem functions, and also that several of the current theoretical ecological hypotheses could coexist to explain the observed soil organic carbon distribution patterns.

CONCLUSION
Our DIM revealed that SOC storage in the Pyrenees not only depends on regional, landscape and local scale factors, including climate and topography, but also on the contribution of the different plant functional guilds and interactions between the guilds. Additionally grass proportion effects were modulated by grazing intensity. In particular, legumes had a complex effect as they enhanced SOC stocks at low proportions and reduced them at high proportions. The magnitude of these effects depended on the relative composition of other guilds in the grassland, being maximum when grass and forb proportions, besides legumes, were codominant. Different effects of the ability of legumes to fix atmospheric N and their high nutrient acquisition strategy were probably behind this pattern. Our results stress the importance of the keystone role of the N 2 fixation rate of legumes on SOC stocks in natural grasslands and provide a strong argument for species diversity conservation efforts under climate change conditions. In addition, our findings can facilitate the elaboration of regional and local strategies to ameliorate the soil capacity to absorb carbon, contributing to the global effort to preserve terrestrial carbon pools.

AC KNOWLEDGEMENT S
We would like to express our thanks to the many people who collaborated in fieldwork, sample processing and data analysis over time. Research in this paper is based on the PASTUS database, compiled from different funding sources over time, the most relevant being: the EU Interreg III-A Pro-gramme (I3A-4-147-E) and the POCTEFA Programme/Interreg IV-A (FLUXPYR, EFA 34/08); the Spanish Science Foundation FECYT-MICINN (CARBOPAS: REN2002-04300-C02-01; CAR-BOAGROPAS: CGL2006-13555-C03-03 and CAPAS: CGL2010-22378-C03-01); the Foundation Catalunya-La Pedrera; and the Spanish Institute of Agronomical Research INIA (CARBOCLUS: SUM2006-00029-C02-0). This work was funded by the Spanish Science Foundation FECYT-MINECO (BIOGEI: GL2013-49142-C2-1-R and IMAGINE: CGL2017-85490-R) and the University of Lleida (PhD Fellowship to AR). This research article has received a grant for its linguistic revision from the Language Institute of the University of Lleida (2021 call).

FUNDING
Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.

OPEN ACCESS
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit h ttp://creativecommons.org/licenses/by/4.0/.