Interactions between biogeochemical and management factors explain soil organic carbon in Pyrenean grasslands

Grasslands are one of the major sinks of terrestrial soil organic carbon (SOC). Understanding how environmental and management factors drive SOC is challenging because they are scale-dependent, with large-scale drivers affecting SOC both directly and through drivers working at small scales. Here we addressed how regional, landscape and grazing management, soil properties and nutrients, and herbage quality factors affect 20 cm depth SOC stocks in mountain grasslands in the Pyrenees. Taking advantage of the high variety of environmental heterogeneity in the Pyrenees, we built a dataset (n=128) that comprises a wide range of environmental and management conditions. This was used to understand the relationship between SOC stocks and their drivers considering multiple environments. We found that temperature seasonality (difference between mean summer temperature and mean annual temperature; TSIS) was the most important geophysical driver of SOC in our study, depending on topography and management. TSIS effects on SOC increased in exposed hillsides, slopy areas, and relatively intensively grazed grasslands. Increased TSIS probably favours plant biomass production, particularly at high altitudes, but landscape and grazing management factors regulate the accumulation of this biomass into SOC. Concerning biochemical SOC drivers, we found unexpected interactive effects between grazer type, soil nutrients and herbage quality. Soil N was a crucial SOC driver as expected but modulated by livestock species and neutral detergent fibre contenting plant biomass; herbage recalcitrance effects varied depending on grazer species. These results highlight the gaps in knowledge about SOC drivers in grasslands under different environmental and management conditions. They may also serve to generate testable hypotheses in later/future studies directed to climate change mitigation policies.

1 Introduction

Soil organic carbon (SOC) is crucial for the functioning of terrestrial ecosystems (Lal, 2004a). SOC enhances soil and water quality and biomass productivity and has an important role in relation to climate change (Lal, 2004b). Although grasslands have small aboveground biomass compared to other ecosystems, their SOC stocks can be comparable to those in forest ecosystems (Berninger et al., 2015). This is due to their high root biomass and residues, which are a substantial carbon source and can contribute to water retention in soil. This creates favourable conditions for the accumulation of organic matter (Von Haden and Dornbush, 2014; Yang et al., 2018). These attributes, together with the high extent of grassland global cover, make grasslands store around 34 % of the terrestrial carbon, mostly in their soils (White et al., 2000).

SOC accumulation results from a complex equilibrium between primary production and organic matter decomposition which depends on multiple environmental factors such as climate, soil texture and nutrients, or land management (Jenny, 1941; Schlesinger, 1977). Understanding how these scale-dependent environmental factors drive SOC is challenging because large-scale drivers also affect those working at fine spatial scales. This has been described as a hierarchy of controls over SOC (Fig. 1; Manning et al., 2015).

Figure 1 Conceptual scheme used in this work to investigate potential environmental drivers with soil organic carbon (SOC). We assume that drivers may affect SOC both directly or hierarchically through another driver. Interactions between factors acting at different scales and belonging to different categories could also drive SOC. Grazing management has a special status because it may be acting at different scales, landscape and local.

Climate is known to be the main SOC driver at broad (global and regional) scales, mean annual precipitation (MAP) and mean temperature (MAT) being the most frequent climate indicators (Wiesmeier et al., 2019). However, climate annual variations represented by seasonality variables are commonly neglected when considering possible SOC drivers in broad-scale models, in spite of being important drivers of plant primary production and enzymatic activity of soil microorganisms (Fernández-Alonso et al., 2018; Garcia-Pausas et al., 2007; Puissant et al., 2018). Other regional and landscape factors like bedrock or topography are also considered to be at the top of the hierarchy because they influence multiple geophysical and biochemical factors affecting SOC, including soil texture and water flow paths (Gray et al., 2015; Hobley et al., 2015). Next in the hierarchy after regional and landscape factors are several soil geophysical properties, like pH and texture, which are controlled by climate and bedrock and which affect SOC through both plant primary production and microbial activity and the capacity to stabilize the SOC (Deng et al., 2016; Xu et al., 2016a).

Soil macro- and micro-nutrients are in the next level of the hierarchy, as their abundance is determined by multiple factors, including climate, soil pH, water content or clay content (Hook and Burke, 2000; de Vries et al., 2012). They play an essential role in primary production and herbage quality and act as resources for microbes to mineralize SOC (Aerts and Chapin, 1999; Vitousek and Howarth, 1991). However, these variables are commonly omitted as possible drivers of SOC in the broad-scale studies, especially in those studies focusing on predictive rather than explicative models (Gray et al., 2015; Manning et al., 2015; Zhang et al., 2018). This kind of variable is less frequently available and more difficult to measure than the other indicators used for SOC modelling (Manning et al., 2015). Furthermore, the use of soil nutrients as SOC drivers in linear models can be challenging, as they are often strongly linked to SOC dynamics. This may mask the effect of other drivers acting at larger spatial scales (Bing et al., 2016; Cleveland and Liptzin, 2007; Tipping et al., 2016).

Vegetation represents another group of SOC drivers affected by climate, topography and soil properties and nutrients (Fernández-Martínez et al., 2014; de Vries et al., 2012; Zhu et al., 2019). Plant biomass is the main input of organic carbon into the soil (Shipley and Parent, 1991). However, a not so frequently considered factor is plant litter quality, which can determine decomposition rates and patterns and hence soil carbon sequestration (Ottoy et al., 2017; Yan et al., 2018, 2019).

In addition to these factors, livestock management effects on grassland SOC are a noteworthy issue since they are poorly understood (Wiesmeier et al., 2019). It is known that herbivores can 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 et al., 2017). Several studies confirmed the interaction between grazing and other SOC drivers at diverse scales (Abdalla et al., 2018; Eze et al., 2018; Lu et al., 2015, 2017; Zhou et al., 2017). Hence, grazing management may be considered a SOC driver with effects at multiple levels of the driver hierarchy (Fig. 1), both affecting other SOC drivers and interacting with them. However, most of the studies investigating grazing effects on SOC focus on grazing intensity in spite of evidence pointing to a greater role of grazer species in determining vegetation and SOC (Chang et al., 2018; Sebastià et al., 2008b).

In this study, our goal was to identify the main drivers of SOC stocks and their interactions in Pyrenean mountain grasslands. For this purpose, we considered a wide set of regional, landscape, soil geophysical and biochemical, and herbage quality factors, together with grazing management factors. Mountain grasslands comprise a wide range of all these conditions, which make carbon stocks highly variable (Garcia-Pausas et al., 2007, 2017). For this reason, data analysed here include a wide range of environmental conditions, comparable to studies on SOC developed at continental or even worldwide scales (Table 1). Additionally, we considered an exceptionally broad compilation of drivers (Table 1). To deal with correlations and interactions between SOC drivers, we developed an exhaustive modelling approach based on the controls over function hypothesis (de Vries et al., 2012). To facilitate the formulation of our specific questions to answer in this study, we classified SOC drivers into three main groups (Fig. 1): (i) geophysical factors, which include regional and landscape factors and are supposed to be the first sources of variation, (ii) biochemical factors, which include soil nutrients and herbage factors and could be conditioned by geophysical factors, and (iii) grazing management factors, which could affect SOC through multiple interactions with the rest of the variables at multiple scales. In particular, the specific questions of this study are the following. (1) What are the relative and interaction effects of the geophysical and biochemical SOC controls? (2) How does grazing management regulate the effects of other SOC drivers?

Table 1 Considered factors affecting SOC stocks in some recent studies. V: the study considers this variable type; –: the study does not consider this variable type. In italics: present study.

¹ It considers SOC concentrations. ² It considers total carbon stocks. ³ It considers total carbon stocks and its fractions. ^a Fertilizer effects. ^b Only aboveground and/or belowground biomass index.

2 Material and methods

2.1 Location and sampling design

The set of data used in this study has been extracted from the PASTUS Database (http://ecofun.ctfc.cat/?p=3538, last access: 12 January 2019), which was 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 128 grassland patches distributed across the Pyrenees (Fig. S1 in the Supplement) and including topographical, climate, soil, herbage and management variables. The elaboration of the PASTUS Database concerning this study is summarized in Fig. 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; Rodríguez et al., 2018; Table 1). Almost all of the plant species in the grasslands from the PASTUS Database are perennial (Sebastià, 2004), and plant diversity is highly heterogeneous, as are the environmental conditions (Rodríguez et al., 2018).

Sampling in the PASTUS Database was designed according to a stratified random scheme, where samples were selected at random within strata. This process was done 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 north-eastern sector and the land use map of Navarra 1 : 25 000 (Gobierno de Navarra, 2003) for the western sectors. 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×10 m² 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² plot (see Rodríguez et al., 2018, for additional details about sampling design). Local variables were assessed inside the 100 m² plots. Aboveground biomass was estimated from herbage cut at ground level in four 50×50 cm² quadrats placed in a 2×2 m² subplot inside the 100 m² plot. Herbage from two of the four quadrats was dried and sent to the laboratory for duplicated chemico-bromatological analysis. In addition, in each quadrat, a 20 cm depth soil core was extracted with a 5×5 cm probe after herbage was removed. The soil sample in the probe was separated into two soil layers: 0–10 and 10–20 cm.

2.2 Regional and landscape environmental drivers

In order to investigate the relationship between SOC and potential environmental drivers, 30 independent environmental variables were initially considered (Table S1 in the Supplement). These variables were grouped into five sets: regional, landscape, livestock management, soil nutrient stocks, and herbage variables.

Regional variables included climate variables and bedrock. Climate variables were determined from Worldclim 2.0 (Fick and Hijmans, 2017). We selected mean annual temperature (MAT), mean summer temperature (MST), mean annual precipitation (MAP) and mean summer precipitation (MSP). The difference between mean summer and mean annual temperature emerged as a relevant explanatory factor of soil organic carbon stocks during previous modelling efforts by one of the co-authors (M.-Teresa Sebastià). Later attempts to improve models by substituting this variable with other temperature indices from climatic databases (Fick and Hijmans, 2017) showed that, for the PASTUS Database, this variable provided higher explanatory power than other temperature seasonality indices. Thus, we decided to keep it, and here we name it the Temperature Seasonality Index of Sebastià (TSIS from now on).

Bedrock type was determined in the field and confirmed by the geographical maps mentioned above. Bedrock was categorized into three categories: basic (marls and calcareous rocks), acidic (mostly sandstones and slates) and heterogeneous.

Landscape variables included topography and soil-type variables. Topography variables included slope, aspect, macrotopography and microtopography. Slope and aspect were determined in the field by clinometer and compass respectively. Macrotopography and microtopography were determined visually in the field. Preliminary modelling efforts suggested the reduction of the four macrotopographical positions initially identified in the field into two: mountain top and south-facing slopes were classified as exposed positions and valley bottoms and north-facing slopes as protected macrotopographical positions. Microtopography included three positions: convexities, concavities and smooth areas. Soil-type variables are described in the following.

2.3 Soil physicochemical analysis

To obtain 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 for further physicochemical analysis. Standard physicochemical soil analyses were performed in the upper 0–10 cm soil layer of all grasslands. Some analyses were also performed on samples from the 10–20 cm soil layer, including soil organic carbon and total nitrogen. For those variables, we combined 0–10 and 10–20 cm values to obtain the whole top 20 cm soil layer.

All soil physicochemical analyses were performed on the fine earth, according to standard soil analysis methods. Textural classes were determined by the Bouyoucos method (Bouyoucos, 1936). Soil pH (measured in water), total organic carbon (TOC), total nitrogen (TN), calcium content (Ca), extractable phosphorus (P), magnesium (Mg) and potassium (K) were measured on air-dried samples (Schöning et al., 2013; Solly et al., 2014). Soil carbonates were determined using the Bernard calcimeter. Total carbon and nitrogen (N) contents of the fine earth were determined by an elemental auto-analyser. The organic C fraction was determined by subtracting inorganic C in the carbonates from the total C. Available phosphorus (P) was extracted by the Olsen method (Olsen, 1954). Ca, Mg and K were extracted by ammonium acetate (Simard, 1993) and measured by flame atomic absorption spectroscopy (AAS); David, 1960). Soil organic carbon (SOC) stocks in the upper 20 cm soil layer were then estimated by 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 et al. (2007). Despite recent studies suggesting that fixed mass SOC stocks estimates are preferable to fixed depth methods because they would be more robust to temporal and land use changes in bulk density (Ellert and Bettany 1995), we chose a fixed depth method for measuring SOC stocks. This decision was based on the fact that our work samples came from natural mountain grasslands, where grazing intensity is always low to moderate, and moreover, herbivore presence is seasonal. Therefore, we do not expect important changes in bulk density due to land use. Additionally fixed mass approaches presented the disadvantages of implying more technical difficulties than fixed depth measures, even in the most modern procedures (Haden et al., 2020), and could not deal well with differences in stoniness.

2.4 Herbage chemical and bromatological analysis and NIRS analysis

All four herbage samples per plot were oven-dried at 60 ^∘C to constant weight to determine aboveground biomass and converted into g m⁻². Two out of the four samples were sent to the laboratory for herbage quality analysis. Dried samples were ground to pass a 1 mm stainless steel screen (Cyclotec 1093 Sample mill, Tecator, Sweden) and stored at 4 ^∘C until it was needed for use.

To develop NIRS prediction models, a random subset of 130 samples was used and analysed in duplicate according to the reference methods mentioned further. Procedures described by AOAC were used to determine dry matter (DM) and ash content or mineral matter (MM). Crude protein (CP) was determined by the Kjeldhal procedure (N × 6.25) using a Kjeltec Auto 1030 Analyser (Tecator, Sweden). Samples were analysed sequentially for neutral detergent fibre (NDF), acid-detergent fibre (ADF) and acid-detergent lignin (ADL) in accordance with the method described in Van Soest et al. (1991) using the Ankom 200 Fibre Analyser incubator (Ankom, USA). The fibre analysis was determined on an ash-free basis and without alpha amylase. We calculated two additional herbage quality indexes often used in the bibliography: NDF∕CP and ADL∕NH (Stockmann et al., 2013). For each subsample the C and N contents (CH and NH) were determined by the Dumas dry combustion method using an Elemental Analyzer EA1108 (Carlo Erba, Milan, Italy).

Afterwards, a total of 200 herbage samples were scanned as described below to collect their NIRS spectra. We estimated chemical and bromatological variables according to the equations derived from the previous calibrations on the initial 130 random samples.

NIRS data were recorded from 1100 to 2500 nm using a FOSS NIRSystem 5000 scanning monochromator (Hillerød, Denmark). Separate calibration equations were generated for grassland samples. Reflectance (R) data were collected in duplicate every 2 nm. WinISI III (v. 1.6) software (FOSS, Denmark) was employed for spectra data analysis and development of chemometric models. Prior to calibration, log 1∕R spectra were corrected for the effects of scatter using the standard normal variate (SNV), detrend (DT) and multiple scatter correction (MSC) and transformed into first or second derivatives using a different gap size (nm) and smoothing interval. For each sample, the mean of the spectra from the two lectures was used. Modified partial least squares (MPLS) was the regression method used for calibration development, and cross-validation was undertaken using the standard methodology in the NIRS software program. The performance of the model was determined by the following statistical tools: standard error of calibration (SEC); standard error of cross-validation (SECV); coefficient of determination for calibration (R²) and cross-validation ($r_{\mathrm{cv}}^{\mathrm{2}}$); the ratio of performance to deviation (RPD) described as the ratio of the standard deviation for the validation samples to the standard error of cross-validation (RPD = SD ∕ SECV) should ideally be at least 3, and the range error ratio (RER = Range ∕ SECV) described as the ratio of the range in the reference data to the SECV should be at least 10 (Williams and Sobering, 1996; Williams et al., 2014).

2.5 Livestock management variables

The management variables (grazer type) initially used for sampling stratification were determined from records available in the municipalities of the study area. Once the specific grassland patches to be sampled were determined, we carried out a detailed analysis of the management where the patches were located. To this effect, we carried out detailed surveys among farmers, shepherds and land managers. Sometimes the information collected was modified according to visual records in the field (e.g. cattle and/or cattle dung found in supposedly ungrazed areas). Information from municipalities was usually the most imprecise.

We considered two management variables: grazing intensity and grazer type. Grazing intensity was determined by estimating livestock stocking rates measured as livestock units ha⁻¹ (LU ha⁻¹) and treated as a semi-quantitative variable with three categories: low (1; lower than 0.2 LU ha⁻¹), medium (2; between 0.2 and 0.4 LU ha⁻¹) and high (3; above 0.4 LU ha⁻¹). Grazer type was categorized into three main types: sheep, cattle and mixed. Mixed grazing included associations comprising small and big livestock species, mainly sheep and cattle, and more rarely horses. Sheep flocks always included some goats.

2.6 Statistical analyses

We applied two different modelling procedures: boosted regression trees (BRTs) and general linear models (GLMs). BRT is an automatic technique that combines insights from traditional statistical modelling and machine learning traditions (Elith et al., 2008). GLMs allowed us to design a hypothesis-based modelling procedure, ensuring that only effects with biological meaning were included, whereas BRT provided information about the data that could be neglected if only a GLM approach was followed.

All the statistical analyses were performed with the R ver. 3.4.3 software (R Core Team, 2016) at a 95 % significance level when appropriate.

2.6.1 Boosted regression trees global model

Including all SOC potential drivers, we fitted a model with BRT to identify the most important variables affecting SOC. BRT uses two algorithms: regression trees and boosting. Regression trees are from the decision tree group of models, and boosting builds and combines a collection of models (Elith et al., 2008). We chose this method because BRT can handle multiple variables better than other techniques such as GLMs and can detect automatically curvilinear relationships and interactions, ignoring non-informative ones. We used the gbm and dismo packages (Greenwell et al., 2019; Hijmans et al., 2017), which provide several functions to fit these models.

Firstly, we fitted a model with all the drivers (Table S1), configured with 15 folds, a Gaussian distribution of the error, a tree complexity of 5, a learning rate of 0.005, a bag fraction of 0.666, and five minimum observations by node. Secondly, we reduced the number of drivers by the method described in Elith et al. (2008). We estimated the change in the model's predictive deviance, dropping one by one each driver, and re-fitted the model with the set of variables which actually improved model performance (Fig. S3). We checked the relative importance of the drivers and the shape and size of the effects by partial effect plots.

2.6.2 General linear models

We designed and executed a modelling procedure based on general linear models (Legendre and Legendre, 1998) and a hierarchy of controls over function (Díaz et al., 2007; de Vries et al., 2012). We log-transformed SOC using the natural logarithm to prevent a breach of the normality assumption by the residuals of the models (Fig. S4). We built two models (Fig. S5), one model based only on geophysical drivers and grazing management (geophysical model) and another model including, in addition to the former drivers, the biochemical drivers: soil nutrients and herbage quality (combined model). With this approach we aimed to avoid ignoring significant effects of the geophysical variables, the original source of variation of SOC stocks according to the hierarchy of controls over function hypothesis (Manning et al., 2015), by masking them with the inclusion of biochemical drivers. We considered that the geophysical factors that potentially affect SOC were regional and landscape (topography and soil-type drivers), as they have been widely used in previous studies to model and predict SOC from landscape to continental scales (Manning et al., 2015; Wiesmeier et al., 2019). In addition to soil nutrients and herbage variables, we included again the livestock management variables in the combined model and looked for interactions involving these variables and biochemical drivers of SOC.

For model building (Fig. S5a), we added driver groups following a sequential order. To fit the geophysical model, we started adding regional, landscape and grazing management drivers and subsequently included soil properties. Afterwards, we sequentially included soil nutrients and herbage drivers to obtain the combined model. We added management variables from the beginning of the modelling process and re-included the discarded ones in each step to guarantee the detection of interactions between management variables and the rest of the drivers. Each time we added a set of drivers, we first considered their main effects and some quadratic terms which were found by preliminary analyses with the scatterplot.matrix function in R package car (Fox et al., 2018); afterwards we included possible level 2 interactions between all the selected drivers.

At every step we selected several candidate terms by a semi-automatic procedure (Fig. S5c) using a genetic algorithm included in R package glmulti (Calcagno, 2015). We used SOC as a response variable in the first step and the residuals of the previous model in the remaining steps (Fig. S5b). This semi-automatic 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 (ΔAICc) values of the best model, as a ΔAICc < 2 indicates that the candidate model is almost as good as the best model (Burnham and Anderson, 2002).

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 modelling procedure (Grueber et al., 2011).

Then, we added these terms to the consolidated model and made a selection through a backward forward procedure. We used several methods to compare and determine the final model, including the AICc, the adjusted determination coefficient R² ($R_{\mathrm{adj}}^{\mathrm{2}}$) and model comparison techniques with the “anova()” function in R, using Chi-square tests to test whether the reduction in the residual sum of squares was statistically significant. Once we had the final model, we assessed the significance of each term by removing it and performing an F test. To estimate 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 (de Vries et al., 2012). We estimated the variance explained by the models through the adjusted determination coefficient R² ($R_{\mathrm{adj}}^{\mathrm{2}}$).

Finally, we estimated the importance of the terms of each model by the lmg method in the relaimpo package (Grömping, 2006) and drew partial effect plots making predictions with R package emmeans (Lenth et al., 2019).

3 Results

SOC stocks of the upper 20 cm layer ranged between 2.6 and 23 kg m⁻², with a median and a mean value of 9.1 and 9.6 kg m⁻² respectively. Standard deviation of the mean was 3.15 (n=125). Minimum, maximum, median and mean values of the continuous drivers are shown in Table S2.

3.1 Relative importance of SOC stock drivers

The final BRT global model achieved a good goodness of fit, with a cross-validated correlation value of 48 % and an explained deviance of 88 %. The most important variables explaining SOC stocks (Fig. 2) were soil N (18 %), soil C∕N (14 %) and clay (13 %), although other variables such as aboveground biomass (7 %), ADL (6 %) or silt (6 %) were also relevant for explaining SOC storage. Three important variables in the BRT model, aboveground biomass, silt and soil K, were not selected in the linear models (Tables 2 and 3). Although accounting for a lower importance value than the previous variables (5 %), TSIS was the most relevant among the climate drivers considered. TSIS was also noticeably important in both linear models (Fig. S6), especially in the geophysical model, not only as the main effect, but also in interaction with other variables (lmg: 4 %–10 %). According to the combined linear model, soil nutrient and herbage variables were other important SOC stock drivers (Fig. S7), but many of these effects occurred in interaction with grazer type.

Figure 2 Relative contributions (%) of driver variables in the final BRT model obtained. Soil N: soil nitrogen; soil C∕N: soil carbon to nitrogen ratio; clay: clay content; abiom: aboveground biomass; ADL: acid-detergent lignin; loam: loam content; K: soil potassium; TSIS: temperature seasonality; NDF: neutro-detergent fibre; pH: soil pH; CH: carbon in the herbage; Mg: soil magnesium; slope: terrain slope; MAP: mean annual precipitation; ADF: acid-detergent fibre. See Table S1 for more details about variables.

Table 2 Results of the geophysical model for soil organic carbon ($R_{\mathrm{adj}}^{\mathrm{2}}=\mathrm{0.34}$). MAP: mean annual precipitation; TSIS: temperature seasonality; slope: terrain slope; exposed: exposed position according to macrotopography; clay: clay content; low and medium intensity: low and medium grazing intensity.

^* P value < 0.05. ${}^{**}$ P value < 0.01. ${}^{***}$ P value < 0.001.

Table 3 Results of the combined model for soil organic carbon ($R_{\mathrm{adj}}^{\mathrm{2}}=\mathrm{0.61}$). MAP: mean annual precipitation; TSIS: mean summer temperature minus mean annual temperature; slope: terrain slope; cattle and mixed: cattle and mixed management according to grazing species; low and medium intensity: low and medium intensity according to grazing intensity; soil C∕N: soil carbon to nitrogen ratio; soil N: soil nitrogen; NDF: neutro-detergent fibre; ADL∕NH: acid-detergent lignin to nitrogen in the herbage ratio.

^* P value < 0.05. ${}^{**}$ P value < 0.01. ${}^{***}$ P value < 0.001.

3.2 Geophysical, biochemical and grazing management effects on SOC stocks

The geophysical model (Table 2) explained 34 % of the total variance ($R_{\mathrm{adj}}^{\mathrm{2}}$). Overall, SOC stocks increased with TSIS under certain conditions: exposed hillsides, high slopes and low stocking rates (Fig. 3a, b and d). On the other hand, clay content had a positive relationship with SOC under low MAP values (Fig. 3c), which became negative at high MAP values.

Figure 3 Relationship between SOC and regional- and landscape-scale factors in the geophysical model. In (a) solid lines and circles represent exposed hillsides, and dotted lines and triangles indicate protected hillsides. In (d) solid lines and circles indicate low grazing intensity, dotted lines and triangles indicate medium grazing management intensity and dashed lines and squares indicate high grazing management intensity. In (a)–(d) line and plane values are predictions of the model across the corresponding predictors' range according to estimated marginal means. Grey areas around regression lines indicate standard errors. In (a) and (d) points indicate actual values.

Adding nutrient and herbage variables to the previous geophysical model to build the combined model (Table 3) increased the total variance ($R_{\mathrm{adj}}^{\mathrm{2}}$) up to 61 %. Macrotopography and clay effects described by the geophysical model were removed by the new model terms included. SOC increased with C∕N (Fig. 4a). Soil nitrogen modulated the effects of livestock type and NDF on SOC. Cattle-grazed grasslands stored more SOC than mixed and sheep-grazed grasslands under low soil N conditions, whereas the reverse occurred at high soil N levels (Fig. 3b). NDF had negative effects on SOC stocks at high soil N values but had no effect under low soil N levels (Fig. 4c). Finally, herbage ADL∕NH had positive effects on SOC under mixed and sheep-grazing regimes, but there was no effect under cattle management (Fig. 4d).

Figure 4 The relationship between SOC and biochemical and herbage factors in the combined model. In (b) and (d) solid lines and circle points represent cattle grazing, dashed lines and square points indicate sheep grazing, and dotted lines and triangle points indicate mixed grazing. In (a)–(d) line and plane values are predictions of the model across the corresponding predictors' range according to estimated marginal means. In (a)–(d) line and plane values are predictions of the model across the corresponding predictors' range according to estimated marginal means. The grey spectrum indicates 95 % confidence intervals. In (a) and (d) points indicate actual values.

4 Discussion

4.1 Considerations about the modelling procedure

Unsurprisingly, the SOC drivers selected and their main effects in both of the modelling approaches (BRTs and GLMs) were highly congruent (Figs. 2–4, S8). Consequently, we preferred to focus on the results from the linear models because this approximation allowed us to build models under a hierarchy of controls over function hypothesis (Manning et al., 2015). Hence, although it is not possible to unequivocally establish the causal links between SOC drivers (Grace, 2006; Grace and Bollen, 2005), with our GLM procedure we guarantee that the effects of the biochemical variables added in the combined model on SOC stocks have not been exclusively induced by geophysical drivers (de Vries et al., 2012). If this was the case, soil nutrient and herbage quality drivers would not have entered the combined model as significant terms. This happened with aboveground biomass, which is assumed to be a very important SOC driver, and indeed aboveground biomass was relevant in the BRT model, but in the GLM it was substituted by other, more meaningful, variables. In addition, our GLM modelling approach enabled us the selection of biologically meaningful interactions (Manning et al., 2015; de Vries et al., 2012), which cannot be done with a fully automatic approach like BRT. This GLM-sequenced modelling procedure, looking for the primary sources of variation, together with the stratified sampling design, is useful as it led us to select a set of lowly correlated drivers for our linear models (Table S5). Furthermore, the BRT model provided some valuable information, identifying some relevant SOC drivers which were discarded during the GML modelling, like aboveground biomass or soil silt and K (Figs. 2 and S8). The effects of those drivers were probably masked by the effects of other variables in our linear models (Yang et al., 2009), indicating that these factors were presumably pathways through which other variables drove SOC (de Vries et al., 2012). These variables, identified by BRT and discarded by GLMs, should be considered as potential SOC drivers in further studies, particularly when more detailed and difficult-to-obtain biochemical variables, present in our database, are not available.

4.2 Geophysical, biochemical and grazing management factors driving SOC stocks

Considering the difficulties of modelling SOC in a widely heterogeneous mountain environment (Garcia-Pausas et al., 2017), the geophysical model provided important information about broad-scale and topographic SOC drivers in the Pyrenees. This information could be useful not only for a better understanding of SOC patterns in mountain grasslands, but also for future modelling studies aiming to predict SOC, since geophysical variables are easier and less expensive to acquire and measure compared to biochemical variables (Manning et al., 2015).

Most studies on soil carbon usually pinpoint mean temperature and precipitation as the most important climate drivers of SOC (Hobley et al., 2015; Manning et al., 2015; Wiesmeier et al., 2019). Climate regulates large-scale patterns of aboveground net primary production (Chapin et al., 1987). In our study, temperature seasonality (TSIS) was a key driver of SOC, modulated by macrotopography, slope and grazing intensity (Table 2; Fig. 3). The highest variation of TSIS in our database, that is, the broadest temperature seasonality, occurred in cold environments, as compared to mild climates (Fig. S9). In mountain grasslands, cold climates imply a short phenological period of development for plants (Gómez, 2008). Hence, the positive effect of TSIS on SOC could be associated with a higher biomass accumulation in cold locations with more favourable temperatures during summer, this fact reducing geophysical stress for plants and broadening their growth period (Garcia-Pausas et al., 2007; Kikvidze et al., 2005). This increase in soil organic matter inputs during summer would overcome an eventual increase in soil organic matter decomposition caused by high temperatures (Sanderman et al., 2003).

The interactive effects of TSIS on SOC stocks with macrotopography and slope illustrate the capacity of landscape factors to modulate macroclimate effects on soil (Hook and Burke, 2000). Induced microclimate changes are often the explanation for the effects of topography in SOC (Lozano-García et al., 2016). In our case, SOC stocks increased with temperature seasonality, particularly in exposed locations, including south-facing hillsides and hillside tops (Fig. 3a; Table 2). In protected locations, including shady hillsides and valley bottoms, the hypothesized positive effect of increased TSIS values on plant productivity could be mitigated due to reduced solar radiation, long snow-covered periods and freezing episodes (Garcia-Pausas et al., 2007; López-Moreno et al., 2013). Additionally, differences in SOC between exposed and protected sites may also occur through other mechanisms, for instance the alteration of soil physicochemical properties (Zhang et al., 2018) or differences in vegetation (Sebastià, 2004). Since we used a hierarchy of controls approach (Manning et al., 2015), these indirect topographical effects on SOC stocks could be behind the exclusion in the linear models of some drivers selected in the BRT model, like silt or pH (Figs. 2 and 3). In addition, SOC stocks decreased with increase in slope, which may be attributed to reduced carbon inputs and increased carbon losses induced by steeper slopes (Yuan et al., 2019, and references therein). However, we found that increased temperature seasonality (TSIS) values partly compensated for negative slope effects on SOC.

The effect of temperature seasonality on SOC stocks was also modified by grazing management. At low TSIS values, SOC stocks increased under moderate to high grazing pressure; this effect disappeared as TSIS values increased (Fig. 3d). Recent meta-analyses concluded that intensive grazing commonly has decreasing effects on SOC (Abdalla et al., 2018; Eze et al., 2018; Mcsherry and Ritchie, 2013). However, these effects were strongly context-specific, depending on other factors, including climate and soil-type vegetation (Abdalla et al., 2018; Eze et al., 2018; Mcsherry and Ritchie, 2013). Moreover, moderate grazing intensities can increase SOC inputs by dung deposition and aboveground and root biomass production (Franzluebbers et al., 2000; Zeng et al., 2015). In our study, grazing intensity was relatively moderate (see methods); therefore, in our study increasing stocking rates may increase soil carbon inputs in moderate seasonality locations by enhancing aboveground and belowground productivity.

Soil texture also showed interactive effects on SOC stocks with climatic variables. In particular, clay effects on SOC stocks became negative as MAP values increased (Fig. 3c; Table 2). Both MAP and clay content are widely assumed to be positively correlated with SOC (Wiesmeier et al., 2019), but high soil water content caused by high MAP may inhibit decomposition if a shortage of oxygen supply occurs (Xu et al., 2016b). Furthermore, fine texture soils could be waterlogged frequently, leading to inhibition of root growth and soil C allocation belowground (Mcsherry and Ritchie, 2013).

4.3 Geophysical, biochemical and grazing management factors driving SOC stocks

Considering the difficulties of modelling SOC in a widely heterogeneous mountain environment (Garcia-Pausas et al., 2017), the geophysical model provided important information about SOC drivers in the Pyrenees. This information could be useful not only for a better understanding of SOC patterns in mountain grasslands, but also for future modelling studies aiming to predict SOC, since geophysical variables are easier and less expensive to acquire and measure compared to biochemical ones (Manning et al., 2015).

TSIS was a key driver of SOC with a varying effect depending on macrotopography, slope and grazing intensity (Table 2; Fig. 3). While most of the previous studies addressing soil carbon did not include any temperature seasonality variable as a potential SOC predictor, usually focusing on mean temperature and precipitation as the most important climate drivers of SOC (Hobley et al., 2015; Manning et al., 2015; Wiesmeier et al., 2019), our models suggest that TSIS and other temperature seasonality indexes should be included in further studies to provide more evidence of the extent of the effects of temperature seasonality on SOC stocks.

Climate regulates large-scale patterns of aboveground net primary production (Chapin et al., 1987). In the case of mountain grasslands, cold climates imply a short phenological period of development for plants (Gómez, 2008). Cold sites characterized by low mean temperatures presented a wider spectrum of TSIS values than warm sites, presenting both the lowest and highest TSIS values (Fig. S9). Hence, the positive effect of TSIS on SOC could be associated with a higher biomass accumulation in cold locations with more favourable temperatures during summer, this fact reducing geophysical stress for plants and broadening their growth period (Garcia-Pausas et al., 2007; Kikvidze et al., 2005). This rise in soil organic matter inputs during summer would overcome an eventual increase in soil organic matter decomposition rates due to high temperatures (Sanderman et al., 2003), which could even be diminished if microbial biomass decreases as a result of soil moisture reduction (Puissant et al., 2018).

The interactions of TSIS with macrotopography and slope illustrate the capacity of landscape factors to modulate macroclimate effects on soil (Hook and Burke, 2000). Induced microclimate changes are often the explanation for the effects of topography in SOC (Lozano-García et al., 2016). In our case, SOC stocks increased with temperature seasonality, particularly at mountain-exposed areas (Fig. 3a; Table 2). In protected sites, located on shady slopes and in valley bottoms, the hypothesized positive effect of high TSIS values on plant productivity could be mitigated due to lower solar radiation, longer snow-covered periods and freezing episodes (Garcia-Pausas et al., 2007; López-Moreno et al., 2013). Conversely, negative effects of low TSIS values on plant productivity could be compensated for thanks to the more humid conditions in protected sites compared to the exposed sites (Garcia-Pausas et al., 2007). Additionally, it is important to take into account that differences in SOC between exposed and protected sites may also occur through other mechanisms, for instance the alteration of soil physicochemical properties like pH, soil texture or stoniness (Zhang et al., 2018) or differences in vegetation (Sebastià, 2004). Since we used a hierarchy of controls approach (Manning et al., 2015), these topography indirect effects on SOC stocks could be behind the exclusion in the linear models of some drivers selected in the BRT model, like silt or pH (Figs. 2 and 3).

In addition, high TSIS values compensated for SOC stock decrease with a greater slope, which may be attributed to reduced carbon inputs and increased carbon losses induced by steeper slopes (Yuan et al., 2019, and references therein). Increases in grazing pressure elevated SOC stocks under low TSIS values (Fig. 3d). This was a surprising result according to recent meta-analyses, which concluded that grazing has commonly decreasing effects on SOC (Abdalla et al., 2018; Eze et al., 2018; Mcsherry and Ritchie, 2013). However, these effects were strongly context-specific, depending on other factors like climate and soil-type vegetation (Abdalla et al., 2018; Eze et al., 2018; Mcsherry and Ritchie, 2013). Moreover, light and medium grazing intensities can increase SOC inputs by dung deposition and promoting aboveground and root biomass production (Franzluebbers et al., 2000; Zeng et al., 2015). Considering that in our natural grasslands all grazing intensities are relatively low (see methods), our medium and high stock rates may increase soil carbon inputs in low-seasonality locations by enhancing aboveground and belowground productivity.

Interestingly, clay content and precipitation presented interacting effects on SOC (Fig. 3c; Table 2). Both MAP and clay content are widely assumed to be positively correlated with SOC (Wiesmeier et al., 2019). High MAP would increase SOC inputs by promoting plant productivity (Jobbágy and Jackson, 2000; Hobley et al., 2015). Clay positive effects are often attributed to a larger contact surface of soil particles (Kennedy et al., 2002), the absorption of negatively charged organic matter, high soil water retention and the exclusion of decomposer organisms due to their low pore size (Krull et al., 2001). In our study, high soil water contents caused by high MAP may inhibit decomposition if a shortage of oxygen supply occurs (Xu et al., 2016b). However, as MAP values increased, the clay effect on SOC became negative. To explain low SOC values at high MAP and high clay content, Mcsherry and Ritchie (2013) hypothesized that finer texture soils could be waterlogged more frequently, leading to inhibition of root growth and soil C allocation belowground.

The addition of soil nutrient and herbage variables to our geophysical model implied substitution of terms, including clay content and macrotopography, by newly added variables (Tables 2 and 3). This highlights the importance of indirect effects of these variables on SOC through other small-scale drivers (Leifeld et al., 2015; Xu et al., 2016b; Zhu et al., 2019). The combined model was complex and included infrequently tested effects involving interactions between grazer type, soil nutrients and herbage quality variables (Table 3, Fig. 4). Those results must be interpreted cautiously, because they are based on observational data but can contribute to generating testable hypotheses for later studies about some complex and untested relationships between SOC and its drivers. Interaction experiments concerning soil properties are expensive and rare in the literature (Rillig et al., 2019).

For this reason, SOC increased with the C∕N ratio (Fig. 4a), which may be explained by the difficulty of soil organic matter decomposition by soil microbes, decreasing decomposition rates of SOC with increasing soil C∕N (Wanyama et al., 2019; Xu et al., 2016b). A positive relationship between SOC and soil N was also expected, since most of the soil N is in combined form with organic matter (Cambardella and Elliott, 1994). However, in this case, due to the wide range of conditions and the randomized sampling design of the PASTUS Database, the raw correlation between soil N and SOC was somehow discrete (r=0.297; p value = 0.001; R²=0.088) in comparison to other studies (i.e. Yan et al., 2020). However, the novelty revealed by our model is that soil N could modulate the effects of certain SOC drivers, including livestock type and herbage NDF.

Cattle-grazed grasslands stored more SOC than mixed and sheep-grazed grasslands, but only under low soil N conditions (Fig. 4b). Chang et al. (2018) found that in a N poor semi-arid grassland, sheep decreased SOC content in comparison to cattle due to vegetation changes caused by their feeding preference for highly palatable forbs (Sebastià et al., 2008a), thus promoting less palatable grasses which supported less root biomass. In overall, under low soil N conditions, palatable plants are expected to contribute to SOC inputs through the stimulation of C allocation in forb roots (Ågren and Franklin, 2003; Warembourg et al., 2003) and the increase in the overall plant productivity due to legume atmospheric N fixation (Van Der Heijden et al., 2008).

However, these processes could decline under high soil N contents. For instance, legume atmospheric N fixation could be reduced since it requires additional energy in comparison to nitrogen acquisition from the soil (Ibañez et al., 2020; Minchin and Witty, 2005). Additionally, sheep selective feeding habits could shift plant leaf traits in the community towards nutrient-conservative leaf traits, which commonly induce fungal-based soil food webs with slow nutrient cycling and high SOC storage due to low decomposition rates (Orwin et al., 2010).

Additionally, grasslands with mixed grazed regimes stored even more SOC than sheep-grazed grasslands under high soil N conditions (Fig. 4b, Table 3). This result emphasizes that mixed livestock assemblages deserve particular attention, because mixed grazing can affect plant composition distinctly from single grazing species regimes and alter traveling and trampling behaviour of grazing animals (Aldezabal et al., 2019; Chang et al., 2018; Liu et al., 2015).

NDF was negatively related to SOC at high soil N values (Fig. 4c). NDF proportion represents the amount of structural compounds on litter and hence is inversely related to non-structural compound content (Goering and Van Soest, 1970). The latter are the main source of organic matter formation at the early stages of decomposition, and they are incorporated into microbial biomass in a highly efficient way (Cotrufo et al., 2013). However, if microbial necromass was recycled by microbes before its incorporation into mineral-associated organic matter (Córdova et al., 2018), it could be respired and mineralized instead of stored. Thus, our model suggests that incorporation of labile and fast metabolized non-organic compounds to soil organic matter could be a pathway of SOC allocation at high soil N in Pyrenean grasslands.

On the other hand, the ADL∕NH ratio was positively related to SOC in sheep and mixed grazed grasslands (Fig. 4d). The ADL∕NH ratio is a commonly used indicator for the resistance of litter to degradation, particularly at later stages of decomposition (Taylor et al., 1989). Hence, the increase in SOC stocks with ADL∕NH should be related to the physical pathway of soil organic matter incorporation, forming coarse particulate organic matter (Cotrufo et al., 2015). Moreover, our model suggests that this pathway would be inhibited under cattle grazing, presumably because of their higher digestive efficiency and thus less recalcitrant faeces (Wang et al., 2018) and their less selective diet compared to sheep, as the latter would avoid plants with a high lignin content, promoting recalcitrant litter (Rosenthal et al., 2012; Sebastià et al., 2008a).

Our results concerning interactions between grazer type and herbage quality provide some evidence of grazing effects not only through alterations of plant communities that were reported by previous studies in the region (Canals and Sebastià, 2000; Sebastià et al., 2008a), but also through interactions with them. Although grazing effects were not the most important factors affecting SOC stocks, this is by far the easiest component to manipulate in order to increase or maintain SOC in soils and face climate change (Komac et al., 2014). Considering our results, we suggest conducting more experiments to investigate grazer-type effects on SOC under different soil nutrient conditions and within plant communities with contrasting herbage quality parameters. Grazing management also has other advantages, such as preventing the accumulation of aboveground C and reducing the risk of forest fires (Nunes and Lourenço, 2017).

One key point of our results reinforces the idea that grazer type might be at least as important as grazing intensity in regulating grassland ecosystem dynamics (Tóth et al., 2018) and highlights the need for a more thorough research effort in disentangling not only grazing intensity, but also grazer-type effects on grassland soil organic carbon and nutrient cycling, under different environmental circumstances. The combined model provided some evidence that grazing may affect SOC not only through alterations of plant communities (Canals and Sebastià, 2000; Sebastià et al., 2008a), but also through interactions with them. Although grazing effects were not the most important factors affecting SOC stocks, this is by far the easiest component to manipulate in order to increase or maintain SOC in soils and face climate change (Komac et al., 2014). Despite the need for precise knowledge on the effects of different land uses on ecosystems for climate change mitigation (Lo et al., 2015), studies addressing grazer-type effects on SOC are scarce (i.e. Zhou et al., 2017; Chang et al., 2018). Considering our results, we suggest conducting more experiments which investigate grazer-type effects on SOC under different soil nutrient conditions and within plant communities with contrasting herbage quality parameters.

5 Conclusion

The models presented here show a series of novel broad-scale and local patterns concerning SOC stocks and their geophysical, biochemical and grazing management drivers. Factors driving SOC stocks often interacted in complex ways, within and between spatio-temporal scales. Temperature seasonality (TSIS) was the most critical geophysical factor, affecting SOC through interactions with topographical drivers and grazing intensity. This illustrates how not only climate mean annual conditions should be considered when modelling SOC drivers, but also seasonal patterns. Concerning biochemical factors, we found that the expected positive effect of soil N was modulated by livestock species and herbage NDF, and herbage recalcitrance effects on SOC varied depending on grazer type. Overall, we found a number of interactions highlighting the need to expand knowledge on grassland SOC drivers under different conditions, especially grazing. The latter is the most easily tractable factor affecting SOC. In conclusion, we provided valuable information for further studies dealing with SOC predictions at multiple broad scales and laid out the basis for generating new testable hypotheses for future studies, which may be useful for designing improved policies to palliate climate change.