Biogeosciences, 17, 6033–6050, 2020 https://doi.org/10.5194/bg-17-6033-2020 © Author(s) 2020. This work is distributed under the Creative Commons Attribution 4.0 License. Interactions between biogeochemical and management factors explain soil organic carbon in Pyrenean grasslands Antonio Rodríguez1,2 , Rosa Maria Canals3 , Josefina Plaixats4 , Elena Albanell4 , Haifa Debouk1,2 , Jordi Garcia-Pausas5 , Leticia San Emeterio6,7 , Àngela Ribas8,9 , Juan José Jimenez10 , and M.-Teresa Sebastià1,2 1 Group GAMES, Department of Horticulture, Botany and Landscaping, School of Agrifood and Forestry Science and Engineering, University of Lleida, Lleida, Spain 2 Laboratory of Functional Ecology and Global Change (ECOFUN), Forest Sciences Centre of Catalonia (CTFC), Solsona, Spain 3 Grupo Ecología y Medio Ambiente and ISFood Institute, Universidad Pública de Navarra, Campus Arrosadia, Pamplona, Spain 4 Group of Ruminant Research (G2R), Department of Animal and Food Sciences, Universitat Autònoma de Barcelona, 08193, Bellaterra, Spain 5 Forest Sciences Centre of Catalonia (CTFC), Solsona, Spain 6 Research Institute on Innovation & Sustainable Development in Food Chain (ISFOOD), Universidad Pública de Navarra, 31006 Pamplona, Spain 7 Departamento de Agronomía, Biotecnología y Alimentación, Universidad Pública de Navarra, 31006 Pamplona, Spain 8 Biologia Animal, Biologia Vegetal i Ecologia, Universitat Autònoma de Barcelona, 08193, Bellaterra, Spain 9 Centre for Ecological Research and Forestry Applications (CREAF), 08193, Bellaterra, Spain 10 ARAID/IPE-CSIC, 22700, Huesca, Spain Correspondence: Antonio Rodríguez (antonio.arharh@gmail.com) Received: 21 February 2020 – Discussion started: 17 March 2020 Revised: 14 October 2020 – Accepted: 22 October 2020 – Published: 3 December 2020 Abstract. Grasslands are one of the major sinks of terres- SOC increased in exposed hillsides, slopy areas, and rela- trial soil organic carbon (SOC). Understanding how environ- tively intensively grazed grasslands. Increased TSIS proba- mental and management factors drive SOC is challenging be- bly favours plant biomass production, particularly at high al- cause they are scale-dependent, with large-scale drivers af- titudes, but landscape and grazing management factors regu- fecting SOC both directly and through drivers working at late the accumulation of this biomass into SOC. Concerning small scales. Here we addressed how regional, landscape biochemical SOC drivers, we found unexpected interactive and grazing management, soil properties and nutrients, and effects between grazer type, soil nutrients and herbage qual- herbage quality factors affect 20 cm depth SOC stocks in ity. Soil N was a crucial SOC driver as expected but modu- mountain grasslands in the Pyrenees. Taking advantage of the lated by livestock species and neutral detergent fibre content- high variety of environmental heterogeneity in the Pyrenees, ing plant biomass; herbage recalcitrance effects varied de- we built a dataset (n = 128) that comprises a wide range of pending on grazer species. These results highlight the gaps in environmental and management conditions. This was used to knowledge about SOC drivers in grasslands under different understand the relationship between SOC stocks and their environmental and management conditions. They may also drivers considering multiple environments. We found that serve to generate testable hypotheses in later/future studies temperature seasonality (difference between mean summer directed to climate change mitigation policies. temperature and mean annual temperature; TSIS) was the most important geophysical driver of SOC in our study, de- pending on topography and management. TSIS effects on Published by Copernicus Publications on behalf of the European Geosciences Union. 6034 A. Rodríguez et al.: Soil organic carbon drivers SOC drivers in broad-scale models, in spite of being impor- tant drivers of plant primary production and enzymatic activ- ity of soil microorganisms (Fernández-Alonso et al., 2018; Garcia-Pausas et al., 2007; Puissant et al., 2018). Other re- gional 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 hier- archy after regional and landscape factors are several soil geophysical properties, like pH and texture, which are con- trolled 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). Figure 1. Conceptual scheme used in this work to investigate po- Soil macro- and micro-nutrients are in the next level of the tential environmental drivers with soil organic carbon (SOC). We hierarchy, as their abundance is determined by multiple fac- assume that drivers may affect SOC both directly or hierarchically tors, including climate, soil pH, water content or clay content through another driver. Interactions between factors acting at dif- ferent scales and belonging to different categories could also drive (Hook and Burke, 2000; de Vries et al., 2012). They play an SOC. Grazing management has a special status because it may be essential role in primary production and herbage quality and acting at different scales, landscape and local. 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 1 Introduction in the broad-scale studies, especially in those studies focus- ing on predictive rather than explicative models (Gray et al., Soil organic carbon (SOC) is crucial for the functioning of 2015; Manning et al., 2015; Zhang et al., 2018). This kind terrestrial ecosystems (Lal, 2004a). SOC enhances soil and of variable is less frequently available and more difficult to water quality and biomass productivity and has an impor- measure than the other indicators used for SOC modelling tant role in relation to climate change (Lal, 2004b). Although (Manning et al., 2015). Furthermore, the use of soil nutri- grasslands have small aboveground biomass compared to ents as SOC drivers in linear models can be challenging, as other ecosystems, their SOC stocks can be comparable to they are often strongly linked to SOC dynamics. This may those in forest ecosystems (Berninger et al., 2015). This is mask the effect of other drivers acting at larger spatial scales due to their high root biomass and residues, which are a sub- (Bing et al., 2016; Cleveland and Liptzin, 2007; Tipping et stantial carbon source and can contribute to water retention in al., 2016). soil. This creates favourable conditions for the accumulation Vegetation represents another group of SOC drivers af- of organic matter (Von Haden and Dornbush, 2014; Yang et fected by climate, topography and soil properties and nutri- al., 2018). These attributes, together with the high extent of ents (Fernández-Martínez et al., 2014; de Vries et al., 2012; grassland global cover, make grasslands store around 34 % Zhu et al., 2019). Plant biomass is the main input of organic of the terrestrial carbon, mostly in their soils (White et al., carbon into the soil (Shipley and Parent, 1991). However, 2000). a not so frequently considered factor is plant litter quality, SOC accumulation results from a complex equilibrium be- which can determine decomposition rates and patterns and tween primary production and organic matter decomposition hence soil carbon sequestration (Ottoy et al., 2017; Yan et which depends on multiple environmental factors such as cli- al., 2018, 2019). mate, soil texture and nutrients, or land management (Jenny, In addition to these factors, livestock management effects 1941; Schlesinger, 1977). Understanding how these scale- on grassland SOC are a noteworthy issue since they are dependent environmental factors drive SOC is challenging poorly understood (Wiesmeier et al., 2019). It is known that because large-scale drivers also affect those working at fine herbivores can affect SOC through different paths, such as spatial scales. This has been described as a hierarchy of con- regulating the quantity and quality of organic matter returned trols over SOC (Fig. 1; Manning et al., 2015). to soil (Bardgett and Wardle, 2003) or affecting soil respi- Climate is known to be the main SOC driver at broad ration and nutrients by animal trampling or soil microbiota (global and regional) scales, mean annual precipitation alteration (Lu et al., 2017). Several studies confirmed the in- (MAP) and mean temperature (MAT) being the most fre- teraction between grazing and other SOC drivers at diverse quent climate indicators (Wiesmeier et al., 2019). However, scales (Abdalla et al., 2018; Eze et al., 2018; Lu et al., 2015, climate annual variations represented by seasonality vari- 2017; Zhou et al., 2017). Hence, grazing management may ables are commonly neglected when considering possible be considered a SOC driver with effects at multiple levels of Biogeosciences, 17, 6033–6050, 2020 https://doi.org/10.5194/bg-17-6033-2020 A. Rodríguez et al.: Soil organic carbon drivers 6035 the driver hierarchy (Fig. 1), both affecting other SOC drivers Rodríguez et al., 2018; Table 1). Almost all of the plant and interacting with them. However, most of the studies in- species in the grasslands from the PASTUS Database are vestigating grazing effects on SOC focus on grazing inten- perennial (Sebastià, 2004), and plant diversity is highly het- sity in spite of evidence pointing to a greater role of grazer erogeneous, as are the environmental conditions (Rodríguez species in determining vegetation and SOC (Chang et al., et al., 2018). 2018; Sebastià et al., 2008b). Sampling in the PASTUS Database was designed accord- In this study, our goal was to identify the main drivers ing to a stratified random scheme, where samples were se- of SOC stocks and their interactions in Pyrenean mountain lected at random within strata. This process was done using grasslands. For this purpose, we considered a wide set of the ArcMap 10 software (ESRI, Redlands, CA, USA). The regional, landscape, soil geophysical and biochemical, and basis for randomization was the map of habitats of Catalonia herbage quality factors, together with grazing management 1 : 50 000 (Carreras and Diego, 2006) for the eastern and cen- factors. Mountain grasslands comprise a wide range of all tral sectors of the Pyrenees, the map of habitats of Madres- these conditions, which make carbon stocks highly variable Coronat 1 : 10 000 (Penin, 1997) for the north-eastern sec- (Garcia-Pausas et al., 2007, 2017). For this reason, data anal- tor and the land use map of Navarra 1 : 25 000 (Gobierno de ysed here include a wide range of environmental conditions, Navarra, 2003) for the western sectors. Four variables were comparable to studies on SOC developed at continental or initially considered for sampling stratification within each even worldwide scales (Table 1). Additionally, we considered sector: altitude (< 1800; 1800–2300; > 2300 m), slope (0– an exceptionally broad compilation of drivers (Table 1). To 20; 20–30; > 30◦ ), macrotopography (mountain top/south- deal with correlations and interactions between SOC drivers, facing slope; valley bottom/north-facing slope) and grazer we developed an exhaustive modelling approach based on type (sheep; cattle; mixed). Accordingly, we determined a the controls over function hypothesis (de Vries et al., 2012). set of homogeneous grassland patches by crossing the strat- To facilitate the formulation of our specific questions to an- ification variable layers. Grassland patches were then listed swer in this study, we classified SOC drivers into three main by type and arranged within each list randomly to determine groups (Fig. 1): (i) geophysical factors, which include re- sampling priority. At least one to two replicates of each patch gional and landscape factors and are supposed to be the first type defined by the stratification variables were sampled. sources of variation, (ii) biochemical factors, which include In each sampled grassland patch, a 10 × 10 m2 plot was soil nutrients and herbage factors and could be conditioned systematically placed in the middle of each homogeneous by geophysical factors, and (iii) grazing management factors, grassland patch, including a particular plant community. We which could affect SOC through multiple interactions with collected soil and vegetation samples and assessed environ- the rest of the variables at multiple scales. In particular, the mental variables within each 100 m2 plot (see Rodríguez et specific questions of this study are the following. (1) What al., 2018, for additional details about sampling design). Lo- are the relative and interaction effects of the geophysical and cal variables were assessed inside the 100 m2 plots. Above- biochemical SOC controls? (2) How does grazing manage- ground biomass was estimated from herbage cut at ground ment regulate the effects of other SOC drivers? level in four 50×50 cm2 quadrats placed in a 2×2 m2 subplot inside the 100 m2 plot. Herbage from two of the four quadrats was dried and sent to the laboratory for duplicated chemico- 2 Material and methods bromatological analysis. In addition, in each quadrat, a 20 cm depth soil core was extracted with a 5 × 5 cm probe after 2.1 Location and sampling design herbage was removed. The soil sample in the probe was sep- arated into two soil layers: 0–10 and 10–20 cm. The set of data used in this study has been extracted from the PASTUS Database (http://ecofun.ctfc.cat/?p=3538, last 2.2 Regional and landscape environmental drivers access: 12 January 2019), which was compiled by the Lab- oratory of Functional Ecology and Global Change (ECO- In order to investigate the relationship between SOC and po- FUN) of the Forest Sciences Centre of Catalonia (CTFC) tential environmental drivers, 30 independent environmental and the University of Lleida (UdL). We sourced a wealth of variables were initially considered (Table S1 in the Supple- data of 128 grassland patches distributed across the Pyrenees ment). These variables were grouped into five sets: regional, (Fig. S1 in the Supplement) and including topographical, cli- landscape, livestock management, soil nutrient stocks, and mate, soil, herbage and management variables. The elabora- herbage variables. tion of the PASTUS Database concerning this study is sum- Regional variables included climate variables and marized in Fig. S2. The sampled area encompasses a wide bedrock. Climate variables were determined from Worldclim variety of temperate and cold–temperate climates with dif- 2.0 (Fick and Hijmans, 2017). We selected mean annual tem- ferent precipitation conditions, depending on altitude and ge- perature (MAT), mean summer temperature (MST), mean ographical location, from Mediterranean to continental and annual precipitation (MAP) and mean summer precipitation Boreo-Alpine environments (de Lamo and Sebastià, 2006; (MSP). The difference between mean summer and mean https://doi.org/10.5194/bg-17-6033-2020 Biogeosciences, 17, 6033–6050, 2020 6036 A. Rodríguez et al.: Soil organic carbon 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. Article Location Lat Long MAP MAT Topography Grazing Soil Soil Herbage (◦ ) (◦ ) (mm) (◦ C) and bedrock Management properties nutrients Present study Pyrenees 42.14 to −1.22 to 964 to 1.1 to V V V V V 43.3 2.26 1586 9.9 Duarte-guardia Worldwide −51.72 to −163.95 to 65 to −21.2 to V – V – Vb et al. (2019) 80.23 158.25 5115 30 Abdalla et al. Worldwide −45.85 to −114 to 150 to 0 to – V V – V (2018) 51 120.7 1650 21 Eze et al. Worldwide −44 to −121 to 120 to −4.8 to – V V Va Vb (2018) 65 175 2000 26.8 Peri et al. Southern −52 to −73.5 to 139 to 4.2 to V V – – V (2018)1 Patagonia −45 65.5 865 11 Zhang et al. Northern 103.5 to 32.5 to 500 to 8 to V V V – – (2018) China 124.16 42.5 1000 14 Zhao et al. Mongolia 41.95 to 108.28 to 150 to −1.3 to – V V – V (2017) 53.93 116.2 400 2.1 Zhou et al. Worldwide −42.1 to −121 to 200 to 0 to – V – – X (2017)2 52.3 175 600 10 Deng et al. Eastern 28.71 to 120.87 to 940 to 16.86 to V – V – X (2016) China 30.45 122.43 1720 18.57 Gray et al. Eastern −16.7 to −31.8 to 500 to 10 to V X X – V (2015) Australia −43.5 −28.7 2000 30 Lu et al. (2017) Qinghai– 27 to 83 to 37 to −4.04 to – V X – – Tibetan 32 108 718 6.3 Plateau Chang et al. Tibet Not Not 397 to 1.7 to V – – – V (2015)1 reported reported 1910 15.5 Manning et al. England 50.77 to −4.43 to 596 to 6.5 to – V V – V 20153 54.58 0.87 3191 10.9 Mcsherry and Worldwide Not Not 180 to Not – V V – V Ritchie (2013) reported reported 950 reported Garcia-Pausas Pyrenees −7 to 42.5 to 1416 to −0.7 to V – V – – et al. (2007) 2.2 42.9 1904 5 1 It considers SOC concentrations. 2 It considers total carbon stocks. 3 It considers total carbon stocks and its fractions. a Fertilizer effects. b Only aboveground and/or belowground biomass index. annual temperature emerged as a relevant explanatory factor rocks), acidic (mostly sandstones and slates) and heteroge- of soil organic carbon stocks during previous modelling neous. efforts by one of the co-authors (M.-Teresa Sebastià). Later Landscape variables included topography and soil-type attempts to improve models by substituting this variable variables. Topography variables included slope, aspect, with other temperature indices from climatic databases (Fick macrotopography and microtopography. Slope and aspect and Hijmans, 2017) showed that, for the PASTUS Database, were determined in the field by clinometer and compass re- this variable provided higher explanatory power than other spectively. Macrotopography and microtopography were de- temperature seasonality indices. Thus, we decided to keep it, termined visually in the field. Preliminary modelling efforts and here we name it the Temperature Seasonality Index of suggested the reduction of the four macrotopographical po- Sebastià (TSIS from now on). sitions initially identified in the field into two: mountain Bedrock type was determined in the field and confirmed by top and south-facing slopes were classified as exposed po- the geographical maps mentioned above. Bedrock was cat- sitions and valley bottoms and north-facing slopes as pro- egorized into three categories: basic (marls and calcareous tected macrotopographical positions. Microtopography in- Biogeosciences, 17, 6033–6050, 2020 https://doi.org/10.5194/bg-17-6033-2020 A. Rodríguez et al.: Soil organic carbon drivers 6037 cluded three positions: convexities, concavities and smooth converted into g m−2 . Two out of the four samples were sent areas. Soil-type variables are described in the following. to the laboratory for herbage quality analysis. Dried samples were ground to pass a 1 mm stainless steel screen (Cyclotec 2.3 Soil physicochemical analysis 1093 Sample mill, Tecator, Sweden) and stored at 4 ◦ C until it was needed for use. To obtain bulk density, we air-dried and weighed the soil To develop NIRS prediction models, a random subset of samples: we then sieved each sample to 2 mm to separate 130 samples was used and analysed in duplicate according stones and gravels from the fine earth fraction. The fine frac- to the reference methods mentioned further. Procedures de- tion was sent to the laboratory for further physicochemical scribed by AOAC were used to determine dry matter (DM) analysis. Standard physicochemical soil analyses were per- and ash content or mineral matter (MM). Crude protein (CP) formed in the upper 0–10 cm soil layer of all grasslands. was determined by the Kjeldhal procedure (N × 6.25) using Some analyses were also performed on samples from the 10– a Kjeltec Auto 1030 Analyser (Tecator, Sweden). Samples 20 cm soil layer, including soil organic carbon and total ni- were analysed sequentially for neutral detergent fibre (NDF), trogen. For those variables, we combined 0–10 and 10–20 cm acid-detergent fibre (ADF) and acid-detergent lignin (ADL) values to obtain the whole top 20 cm soil layer. in accordance with the method described in Van Soest et All soil physicochemical analyses were performed on al. (1991) using the Ankom 200 Fibre Analyser incubator the fine earth, according to standard soil analysis meth- (Ankom, USA). The fibre analysis was determined on an ash- ods. Textural classes were determined by the Bouyoucos free basis and without alpha amylase. We calculated two ad- method (Bouyoucos, 1936). Soil pH (measured in water), to- ditional herbage quality indexes often used in the bibliogra- tal organic carbon (TOC), total nitrogen (TN), calcium con- phy: NDF/CP and ADL/NH (Stockmann et al., 2013). For tent (Ca), extractable phosphorus (P), magnesium (Mg) and each subsample the C and N contents (CH and NH) were potassium (K) were measured on air-dried samples (Schön- determined by the Dumas dry combustion method using an ing et al., 2013; Solly et al., 2014). Soil carbonates were Elemental Analyzer EA1108 (Carlo Erba, Milan, Italy). determined using the Bernard calcimeter. Total carbon and Afterwards, a total of 200 herbage samples were scanned nitrogen (N) contents of the fine earth were determined by as described below to collect their NIRS spectra. We esti- an elemental auto-analyser. The organic C fraction was de- mated chemical and bromatological variables according to termined by subtracting inorganic C in the carbonates from the equations derived from the previous calibrations on the the total C. Available phosphorus (P) was extracted by the initial 130 random samples. Olsen method (Olsen, 1954). Ca, Mg and K were extracted NIRS data were recorded from 1100 to 2500 nm using a by ammonium acetate (Simard, 1993) and measured by flame FOSS NIRSystem 5000 scanning monochromator (Hillerød, atomic absorption spectroscopy (AAS); David, 1960). Soil Denmark). Separate calibration equations were generated for organic carbon (SOC) stocks in the upper 20 cm soil layer grassland samples. Reflectance (R) data were collected in du- were then estimated by taking into account the organic C plicate every 2 nm. WinISI III (v. 1.6) software (FOSS, Den- concentration in the sample and its bulk density and subtract- mark) was employed for spectra data analysis and develop- ing the coarse-particle (> 2 mm) content, following García- ment of chemometric models. Prior to calibration, log 1/R Pausas et al. (2007). Despite recent studies suggesting that spectra were corrected for the effects of scatter using the fixed mass SOC stocks estimates are preferable to fixed depth standard normal variate (SNV), detrend (DT) and multiple methods because they would be more robust to temporal and scatter correction (MSC) and transformed into first or second land use changes in bulk density (Ellert and Bettany 1995), derivatives using a different gap size (nm) and smoothing in- we chose a fixed depth method for measuring SOC stocks. terval. For each sample, the mean of the spectra from the two This decision was based on the fact that our work samples lectures was used. Modified partial least squares (MPLS) was came from natural mountain grasslands, where grazing in- the regression method used for calibration development, and tensity is always low to moderate, and moreover, herbivore cross-validation was undertaken using the standard method- presence is seasonal. Therefore, we do not expect important ology in the NIRS software program. The performance of changes in bulk density due to land use. Additionally fixed the model was determined by the following statistical tools: mass approaches presented the disadvantages of implying standard error of calibration (SEC); standard error of cross- more technical difficulties than fixed depth measures, even in validation (SECV); coefficient of determination for calibra- the most modern procedures (Haden et al., 2020), and could tion (R 2 ) and cross-validation (rcv 2 ); the ratio of performance not deal well with differences in stoniness. to deviation (RPD) described as the ratio of the standard devi- ation for the validation samples to the standard error of cross- 2.4 Herbage chemical and bromatological analysis and validation (RPD = SD / SECV) should ideally be at least 3, NIRS analysis and the range error ratio (RER = Range / SECV) described as the ratio of the range in the reference data to the SECV All four herbage samples per plot were oven-dried at 60 ◦ C should be at least 10 (Williams and Sobering, 1996; Williams to constant weight to determine aboveground biomass and et al., 2014). https://doi.org/10.5194/bg-17-6033-2020 Biogeosciences, 17, 6033–6050, 2020 6038 A. Rodríguez et al.: Soil organic carbon drivers 2.5 Livestock management variables ror, a tree complexity of 5, a learning rate of 0.005, a bag fraction of 0.666, and five minimum observations by node. The management variables (grazer type) initially used for Secondly, we reduced the number of drivers by the method sampling stratification were determined from records avail- described in Elith et al. (2008). We estimated the change in able in the municipalities of the study area. Once the specific the model’s predictive deviance, dropping one by one each grassland patches to be sampled were determined, we carried driver, and re-fitted the model with the set of variables which out a detailed analysis of the management where the patches actually improved model performance (Fig. S3). We checked were located. To this effect, we carried out detailed surveys the relative importance of the drivers and the shape and size among farmers, shepherds and land managers. Sometimes of the effects by partial effect plots. the information collected was modified according to visual records in the field (e.g. cattle and/or cattle dung found in 2.6.2 General linear models supposedly ungrazed areas). Information from municipalities was usually the most imprecise. We designed and executed a modelling procedure based on We considered two management variables: grazing inten- general linear models (Legendre and Legendre, 1998) and sity and grazer type. Grazing intensity was determined by a hierarchy of controls over function (Díaz et al., 2007; de estimating livestock stocking rates measured as livestock Vries et al., 2012). We log-transformed SOC using the natu- units ha−1 (LU ha−1 ) and treated as a semi-quantitative vari- ral logarithm to prevent a breach of the normality assumption able with three categories: low (1; lower than 0.2 LU ha−1 ), by the residuals of the models (Fig. S4). We built two models medium (2; between 0.2 and 0.4 LU ha−1 ) and high (3; above (Fig. S5), one model based only on geophysical drivers and 0.4 LU ha−1 ). Grazer type was categorized into three main grazing management (geophysical model) and another model types: sheep, cattle and mixed. Mixed grazing included asso- including, in addition to the former drivers, the biochemical ciations comprising small and big livestock species, mainly drivers: soil nutrients and herbage quality (combined model). sheep and cattle, and more rarely horses. Sheep flocks always With this approach we aimed to avoid ignoring significant included some goats. effects of the geophysical variables, the original source of variation of SOC stocks according to the hierarchy of con- 2.6 Statistical analyses trols over function hypothesis (Manning et al., 2015), by masking them with the inclusion of biochemical drivers. We We applied two different modelling procedures: boosted re- considered that the geophysical factors that potentially affect gression trees (BRTs) and general linear models (GLMs). SOC were regional and landscape (topography and soil-type BRT is an automatic technique that combines insights from drivers), as they have been widely used in previous studies to traditional statistical modelling and machine learning tra- model and predict SOC from landscape to continental scales ditions (Elith et al., 2008). GLMs allowed us to design a (Manning et al., 2015; Wiesmeier et al., 2019). In addition hypothesis-based modelling procedure, ensuring that only ef- to soil nutrients and herbage variables, we included again fects with biological meaning were included, whereas BRT the livestock management variables in the combined model provided information about the data that could be neglected and looked for interactions involving these variables and bio- if only a GLM approach was followed. chemical drivers of SOC. All the statistical analyses were performed with the R ver. For model building (Fig. S5a), we added driver groups 3.4.3 software (R Core Team, 2016) at a 95 % significance following a sequential order. To fit the geophysical model, level when appropriate. we started adding regional, landscape and grazing manage- ment drivers and subsequently included soil properties. Af- 2.6.1 Boosted regression trees global model terwards, we sequentially included soil nutrients and herbage drivers to obtain the combined model. We added manage- Including all SOC potential drivers, we fitted a model with ment variables from the beginning of the modelling process BRT to identify the most important variables affecting SOC. and re-included the discarded ones in each step to guaran- BRT uses two algorithms: regression trees and boosting. Re- tee the detection of interactions between management vari- gression trees are from the decision tree group of models, and ables and the rest of the drivers. Each time we added a set boosting builds and combines a collection of models (Elith of drivers, we first considered their main effects and some et al., 2008). We chose this method because BRT can handle quadratic terms which were found by preliminary analyses multiple variables better than other techniques such as GLMs with the scatterplot.matrix function in R package car (Fox et and can detect automatically curvilinear relationships and in- al., 2018); afterwards we included possible level 2 interac- teractions, ignoring non-informative ones. We used the gbm tions between all the selected drivers. and dismo packages (Greenwell et al., 2019; Hijmans et al., At every step we selected several candidate terms by a 2017), which provide several functions to fit these models. semi-automatic procedure (Fig. S5c) using a genetic algo- Firstly, we fitted a model with all the drivers (Table S1), rithm included in R package glmulti (Calcagno, 2015). We configured with 15 folds, a Gaussian distribution of the er- used SOC as a response variable in the first step and the resid- Biogeosciences, 17, 6033–6050, 2020 https://doi.org/10.5194/bg-17-6033-2020 A. Rodríguez et al.: Soil organic carbon drivers 6039 uals of the previous model in the remaining steps (Fig. S5b). This semi-automatic process began by obtaining a best sub- set of models using the corrected Akaike information crite- rion (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 al- gorithm, which were those within two Akaike information criterion-corrected (1AICc) values of the best model, as a 1AICc < 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 un- standardized coefficient in the model by the partial standard Figure 2. Relative contributions (%) of driver variables in the final deviation of the variable, which is a function of the standard BRT model obtained. Soil N: soil nitrogen; soil C/N: soil carbon deviation of the variable, the sample size, the number of vari- to nitrogen ratio; clay: clay content; abiom: aboveground biomass; ADL: acid-detergent lignin; loam: loam content; K: soil potassium; ables in the model and the variance inflation factor of the TSIS: temperature seasonality; NDF: neutro-detergent fibre; pH: variable (Barton, 2015). We selected all the variables with soil pH; CH: carbon in the herbage; Mg: soil magnesium; slope: ter- significant effects (alone or in interaction with each other) rain slope; MAP: mean annual precipitation; ADF: acid-detergent in the conditional average model, which was preferred over fibre. See Table S1 for more details about variables. the full average model because we wanted to avoid excessive shrinkage effects at this moment of the modelling procedure (Grueber et al., 2011). 3.1 Relative importance of SOC stock drivers 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 The final BRT global model achieved a good goodness of fit, model, including the AICc, the adjusted determination co- with a cross-validated correlation value of 48 % and an ex- 2 ) and model comparison techniques with efficient R 2 (Radj plained deviance of 88 %. The most important variables ex- the “anova()” function in R, using Chi-square tests to test plaining SOC stocks (Fig. 2) were soil N (18 %), soil C/N whether the reduction in the residual sum of squares was sta- (14 %) and clay (13 %), although other variables such as tistically significant. Once we had the final model, we as- aboveground biomass (7 %), ADL (6 %) or silt (6 %) were sessed the significance of each term by removing it and per- also relevant for explaining SOC storage. Three important forming an F test. To estimate the significance of the main variables in the BRT model, aboveground biomass, silt and effects, we also removed the interaction terms in which they soil K, were not selected in the linear models (Tables 2 and were involved to avoid transferring the effects of the main 3). Although accounting for a lower importance value than terms to the interaction terms (de Vries et al., 2012). We esti- the previous variables (5 %), TSIS was the most relevant mated the variance explained by the models through the ad- among the climate drivers considered. TSIS was also notice- 2 ). justed determination coefficient R 2 (Radj ably important in both linear models (Fig. S6), especially in Finally, we estimated the importance of the terms of each the geophysical model, not only as the main effect, but also model by the lmg method in the relaimpo package (Grömp- in interaction with other variables (lmg: 4 %–10 %). Accord- ing, 2006) and drew partial effect plots making predictions ing to the combined linear model, soil nutrient and herbage with R package emmeans (Lenth et al., 2019). variables were other important SOC stock drivers (Fig. S7), but many of these effects occurred in interaction with grazer type. 3 Results 3.2 Geophysical, biochemical and grazing management SOC stocks of the upper 20 cm layer ranged between 2.6 effects on SOC stocks and 23 kg m−2 , with a median and a mean value of 9.1 and 9.6 kg m−2 respectively. Standard deviation of the mean was The geophysical model (Table 2) explained 34 % of the to- 3.15 (n = 125). Minimum, maximum, median and mean val- 2 ). Overall, SOC stocks increased with TSIS tal variance (Radj ues of the continuous drivers are shown in Table S2. 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. https://doi.org/10.5194/bg-17-6033-2020 Biogeosciences, 17, 6033–6050, 2020 6040 A. Rodríguez et al.: Soil organic carbon drivers 2 = 0.34). MAP: mean annual precipitation; TSIS: temperature sea- Table 2. Results of the geophysical model for soil organic carbon (Radj sonality; slope: terrain slope; exposed: exposed position according to macrotopography; clay: clay content; low and medium intensity: low and medium grazing intensity. Model term Estimate SE t value P value Intercept −0.525 1.802 −0.291 0.771 Climate variables MAP 0.003 0.001 4.560 < 0.001 ∗∗∗ TSIS −0.098 0.228 −0.429 0.669 Topography variables Slope −0.339 0.095 −3.569 0.001 ∗∗∗ Exposed −3.130 0.936 −3.344 0.001 ∗∗ Soil-type variables Clay 0.121 0.027 4.500 < 0.001 ∗∗∗ Management variables Low intensity −5.013 1.196 −4.192 < 0.001 ∗∗∗ Medium intensity 2.012 1.168 1.722 0.088 Interactions TSIS × Exposed 0.417 0.124 3.358 0.001 ∗∗ TSIS × Slope 0.044 0.013 3.452 0.001 ∗∗∗ MAP × Clay 9 × 10−5 2 × 10−5 −4.637 < 0.001 ∗∗∗ TSIS × Low intensity 0.655 0.159 4.110 < 0.001 ∗∗∗ TSIS × Medium intensity −0.262 0.156 −1.684 0.095 ∗ P value < 0.05. ∗∗ P value < 0.01. ∗∗∗ P value < 0.001. Adding nutrient and herbage variables to the previous geo- al., 2015). Hence, although it is not possible to unequivo- physical model to build the combined model (Table 3) in- cally establish the causal links between SOC drivers (Grace, 2 ) up to 61 %. Macrotopog- creased the total variance (Radj 2006; Grace and Bollen, 2005), with our GLM procedure raphy and clay effects described by the geophysical model we guarantee that the effects of the biochemical variables were removed by the new model terms included. SOC in- added in the combined model on SOC stocks have not been creased with C/N (Fig. 4a). Soil nitrogen modulated the ef- exclusively induced by geophysical drivers (de Vries et al., fects of livestock type and NDF on SOC. Cattle-grazed grass- 2012). If this was the case, soil nutrient and herbage qual- lands stored more SOC than mixed and sheep-grazed grass- ity drivers would not have entered the combined model as lands under low soil N conditions, whereas the reverse oc- significant terms. This happened with aboveground biomass, curred at high soil N levels (Fig. 3b). NDF had negative ef- which is assumed to be a very important SOC driver, and in- fects on SOC stocks at high soil N values but had no effect deed aboveground biomass was relevant in the BRT model, under low soil N levels (Fig. 4c). Finally, herbage ADL/NH but in the GLM it was substituted by other, more meaning- had positive effects on SOC under mixed and sheep-grazing ful, variables. In addition, our GLM modelling approach en- regimes, but there was no effect under cattle management abled us the selection of biologically meaningful interactions (Fig. 4d). (Manning et al., 2015; de Vries et al., 2012), which can- not be done with a fully automatic approach like BRT. This GLM-sequenced modelling procedure, looking for the pri- 4 Discussion mary sources of variation, together with the stratified sam- pling design, is useful as it led us to select a set of lowly 4.1 Considerations about the modelling procedure correlated drivers for our linear models (Table S5). Further- more, the BRT model provided some valuable information, Unsurprisingly, the SOC drivers selected and their main ef- identifying some relevant SOC drivers which were discarded fects in both of the modelling approaches (BRTs and GLMs) during the GML modelling, like aboveground biomass or soil were highly congruent (Figs. 2–4, S8). Consequently, we pre- silt and K (Figs. 2 and S8). The effects of those drivers were ferred to focus on the results from the linear models because probably masked by the effects of other variables in our lin- this approximation allowed us to build models under a hi- ear models (Yang et al., 2009), indicating that these factors erarchy of controls over function hypothesis (Manning et Biogeosciences, 17, 6033–6050, 2020 https://doi.org/10.5194/bg-17-6033-2020 A. Rodríguez et al.: Soil organic carbon drivers 6041 Table 3. Results of the combined model for soil organic carbon (Radj 2 = 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. Model term Estimate SE t value P value Intercept −0.290 1.458 −0.199 0.843 Climate variables MAP −0.001 0.000 −2.434 0.017 ∗ TSIS −0.004 0.181 −0.022 0.982 Topography variables Slope −0.225 0.078 −2.868 0.005 ∗∗ Management variables Cattle 0.487 0.101 4.834 < 0.001 ∗∗∗ Mixed −0.289 0.093 −3.106 0.002 ∗∗ Low intensity −3.249 1.014 −3.204 0.002 ∗∗ Medium intensity 1.666 1.073 1.553 0.123 Soil nutrient variables Log(Soil C/N) 0.665 0.076 8.777 < 0.001 ∗∗∗ Soil N 3.302 0.617 5.349 < 0.001 ∗∗∗ Herbage variables NDF 0.014 0.006 2.525 0.013 ∗ Herbage ADL/NH 0.026 0.009 2.987 0.003 ∗∗ Interactions between variable types TSIS × Slope 0.030 0.010 2.833 0.006 ∗∗ TSIS × Low intensity 0.423 0.136 3.104 0.002 ∗∗ TSIS × Medium intensity −0.214 0.143 −1.495 0.138 Soil N × Cattle grazing −0.736 0.168 −4.380 < 0.001 ∗∗∗ Soil N × Mixed grazing 0.493 0.175 2.813 0.006 ∗∗ Soil N × NDF −0.039 0.011 −3.505 0.001 ∗∗∗ Cattle × Herbage ADL/NH −0.030 0.010 −2.872 0.005 ∗∗ Mixed × Herbage ADL/NH 0.014 0.011 1.252 0.213 ∗ P value < 0.05. ∗∗ P value < 0.01. ∗∗∗ P value < 0.001. were presumably pathways through which other variables to acquire and measure compared to biochemical variables drove SOC (de Vries et al., 2012). These variables, identi- (Manning et al., 2015). fied by BRT and discarded by GLMs, should be considered Most studies on soil carbon usually pinpoint mean temper- as potential SOC drivers in further studies, particularly when ature and precipitation as the most important climate drivers more detailed and difficult-to-obtain biochemical variables, of SOC (Hobley et al., 2015; Manning et al., 2015; Wies- present in our database, are not available. meier et al., 2019). Climate regulates large-scale patterns of aboveground net primary production (Chapin et al., 1987). In 4.2 Geophysical, biochemical and grazing management our study, temperature seasonality (TSIS) was a key driver of factors driving SOC stocks SOC, modulated by macrotopography, slope and grazing in- tensity (Table 2; Fig. 3). The highest variation of TSIS in our Considering the difficulties of modelling SOC in a widely database, that is, the broadest temperature seasonality, oc- heterogeneous mountain environment (Garcia-Pausas et al., curred in cold environments, as compared to mild climates 2017), the geophysical model provided important informa- (Fig. S9). In mountain grasslands, cold climates imply a tion about broad-scale and topographic SOC drivers in the short phenological period of development for plants (Gómez, Pyrenees. This information could be useful not only for a bet- 2008). Hence, the positive effect of TSIS on SOC could be ter understanding of SOC patterns in mountain grasslands, associated with a higher biomass accumulation in cold loca- but also for future modelling studies aiming to predict SOC, tions with more favourable temperatures during summer, this since geophysical variables are easier and less expensive fact reducing geophysical stress for plants and broadening https://doi.org/10.5194/bg-17-6033-2020 Biogeosciences, 17, 6033–6050, 2020 6042 A. Rodríguez et al.: Soil organic carbon drivers 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. their growth period (Garcia-Pausas et al., 2007; Kikvidze et controls approach (Manning et al., 2015), these indirect to- al., 2005). This increase in soil organic matter inputs during pographical effects on SOC stocks could be behind the ex- summer would overcome an eventual increase in soil organic clusion in the linear models of some drivers selected in the matter decomposition caused by high temperatures (Sander- BRT model, like silt or pH (Figs. 2 and 3). In addition, SOC man et al., 2003). stocks decreased with increase in slope, which may be at- The interactive effects of TSIS on SOC stocks with macro- tributed to reduced carbon inputs and increased carbon losses topography and slope illustrate the capacity of landscape fac- induced by steeper slopes (Yuan et al., 2019, and references tors to modulate macroclimate effects on soil (Hook and therein). However, we found that increased temperature sea- Burke, 2000). Induced microclimate changes are often the sonality (TSIS) values partly compensated for negative slope explanation for the effects of topography in SOC (Lozano- effects on SOC. García et al., 2016). In our case, SOC stocks increased with The effect of temperature seasonality on SOC stocks was temperature seasonality, particularly in exposed locations, also modified by grazing management. At low TSIS val- including south-facing hillsides and hillside tops (Fig. 3a; ues, SOC stocks increased under moderate to high grazing Table 2). In protected locations, including shady hillsides pressure; this effect disappeared as TSIS values increased and valley bottoms, the hypothesized positive effect of in- (Fig. 3d). Recent meta-analyses concluded that intensive creased TSIS values on plant productivity could be mitigated grazing commonly has decreasing effects on SOC (Abdalla due to reduced solar radiation, long snow-covered periods et al., 2018; Eze et al., 2018; Mcsherry and Ritchie, 2013). and freezing episodes (Garcia-Pausas et al., 2007; López- However, these effects were strongly context-specific, de- Moreno et al., 2013). Additionally, differences in SOC be- pending on other factors, including climate and soil-type veg- tween exposed and protected sites may also occur through etation (Abdalla et al., 2018; Eze et al., 2018; Mcsherry and other mechanisms, for instance the alteration of soil physic- Ritchie, 2013). Moreover, moderate grazing intensities can ochemical properties (Zhang et al., 2018) or differences in increase SOC inputs by dung deposition and aboveground vegetation (Sebastià, 2004). Since we used a hierarchy of and root biomass production (Franzluebbers et al., 2000; Biogeosciences, 17, 6033–6050, 2020 https://doi.org/10.5194/bg-17-6033-2020 A. Rodríguez et al.: Soil organic carbon drivers 6043 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. Zeng et al., 2015). In our study, grazing intensity was rel- could be useful not only for a better understanding of SOC atively moderate (see methods); therefore, in our study in- patterns in mountain grasslands, but also for future modelling creasing stocking rates may increase soil carbon inputs in studies aiming to predict SOC, since geophysical variables moderate seasonality locations by enhancing aboveground are easier and less expensive to acquire and measure com- and belowground productivity. pared to biochemical ones (Manning et al., 2015). Soil texture also showed interactive effects on SOC stocks TSIS was a key driver of SOC with a varying effect de- with climatic variables. In particular, clay effects on SOC pending on macrotopography, slope and grazing intensity stocks became negative as MAP values increased (Fig. 3c; (Table 2; Fig. 3). While most of the previous studies address- Table 2). Both MAP and clay content are widely assumed to ing soil carbon did not include any temperature seasonality be positively correlated with SOC (Wiesmeier et al., 2019), variable as a potential SOC predictor, usually focusing on but high soil water content caused by high MAP may inhibit mean temperature and precipitation as the most important decomposition if a shortage of oxygen supply occurs (Xu et climate drivers of SOC (Hobley et al., 2015; Manning et al., al., 2016b). Furthermore, fine texture soils could be water- 2015; Wiesmeier et al., 2019), our models suggest that TSIS logged frequently, leading to inhibition of root growth and and other temperature seasonality indexes should be included soil C allocation belowground (Mcsherry and Ritchie, 2013). in further studies to provide more evidence of the extent of the effects of temperature seasonality on SOC stocks. 4.3 Geophysical, biochemical and grazing management Climate regulates large-scale patterns of aboveground net factors driving SOC stocks primary production (Chapin et al., 1987). In the case of mountain grasslands, cold climates imply a short phenolog- Considering the difficulties of modelling SOC in a widely ical period of development for plants (Gómez, 2008). Cold heterogeneous mountain environment (Garcia-Pausas et al., sites characterized by low mean temperatures presented a 2017), the geophysical model provided important informa- wider spectrum of TSIS values than warm sites, presenting tion about SOC drivers in the Pyrenees. This information both the lowest and highest TSIS values (Fig. S9). Hence, https://doi.org/10.5194/bg-17-6033-2020 Biogeosciences, 17, 6033–6050, 2020 6044 A. Rodríguez et al.: Soil organic carbon drivers the positive effect of TSIS on SOC could be associated with in low-seasonality locations by enhancing aboveground and a higher biomass accumulation in cold locations with more belowground productivity. favourable temperatures during summer, this fact reducing Interestingly, clay content and precipitation presented in- geophysical stress for plants and broadening their growth pe- teracting effects on SOC (Fig. 3c; Table 2). Both MAP and riod (Garcia-Pausas et al., 2007; Kikvidze et al., 2005). This clay content are widely assumed to be positively correlated rise in soil organic matter inputs during summer would over- with SOC (Wiesmeier et al., 2019). High MAP would in- come an eventual increase in soil organic matter decomposi- crease SOC inputs by promoting plant productivity (Jobbágy tion rates due to high temperatures (Sanderman et al., 2003), and Jackson, 2000; Hobley et al., 2015). Clay positive effects which could even be diminished if microbial biomass de- are often attributed to a larger contact surface of soil particles creases as a result of soil moisture reduction (Puissant et al., (Kennedy et al., 2002), the absorption of negatively charged 2018). organic matter, high soil water retention and the exclusion of The interactions of TSIS with macrotopography and slope decomposer organisms due to their low pore size (Krull et illustrate the capacity of landscape factors to modulate al., 2001). In our study, high soil water contents caused by macroclimate effects on soil (Hook and Burke, 2000). In- high MAP may inhibit decomposition if a shortage of oxy- duced microclimate changes are often the explanation for the gen supply occurs (Xu et al., 2016b). However, as MAP val- effects of topography in SOC (Lozano-García et al., 2016). ues increased, the clay effect on SOC became negative. To In our case, SOC stocks increased with temperature season- explain low SOC values at high MAP and high clay content, ality, particularly at mountain-exposed areas (Fig. 3a; Ta- Mcsherry and Ritchie (2013) hypothesized that finer texture ble 2). In protected sites, located on shady slopes and in val- soils could be waterlogged more frequently, leading to inhi- ley bottoms, the hypothesized positive effect of high TSIS bition of root growth and soil C allocation belowground. values on plant productivity could be mitigated due to lower The addition of soil nutrient and herbage variables to our solar radiation, longer snow-covered periods and freezing geophysical model implied substitution of terms, including episodes (Garcia-Pausas et al., 2007; López-Moreno et al., clay content and macrotopography, by newly added variables 2013). Conversely, negative effects of low TSIS values on (Tables 2 and 3). This highlights the importance of indirect plant productivity could be compensated for thanks to the effects of these variables on SOC through other small-scale more humid conditions in protected sites compared to the drivers (Leifeld et al., 2015; Xu et al., 2016b; Zhu et al., exposed sites (Garcia-Pausas et al., 2007). Additionally, it is 2019). The combined model was complex and included infre- important to take into account that differences in SOC be- quently tested effects involving interactions between grazer tween exposed and protected sites may also occur through type, soil nutrients and herbage quality variables (Table 3, other mechanisms, for instance the alteration of soil physico- Fig. 4). Those results must be interpreted cautiously, because chemical properties like pH, soil texture or stoniness (Zhang they are based on observational data but can contribute to et al., 2018) or differences in vegetation (Sebastià, 2004). generating testable hypotheses for later studies about some Since we used a hierarchy of controls approach (Manning et complex and untested relationships between SOC and its al., 2015), these topography indirect effects on SOC stocks drivers. Interaction experiments concerning soil properties could be behind the exclusion in the linear models of some are expensive and rare in the literature (Rillig et al., 2019). drivers selected in the BRT model, like silt or pH (Figs. 2 and For this reason, SOC increased with the C/N ratio 3). (Fig. 4a), which may be explained by the difficulty of soil In addition, high TSIS values compensated for SOC stock organic matter decomposition by soil microbes, decreas- decrease with a greater slope, which may be attributed to re- ing decomposition rates of SOC with increasing soil C/N duced carbon inputs and increased carbon losses induced by (Wanyama et al., 2019; Xu et al., 2016b). A positive rela- steeper slopes (Yuan et al., 2019, and references therein). In- tionship between SOC and soil N was also expected, since creases in grazing pressure elevated SOC stocks under low most of the soil N is in combined form with organic mat- TSIS values (Fig. 3d). This was a surprising result according ter (Cambardella and Elliott, 1994). However, in this case, to recent meta-analyses, which concluded that grazing has due to the wide range of conditions and the randomized sam- commonly decreasing effects on SOC (Abdalla et al., 2018; pling design of the PASTUS Database, the raw correlation Eze et al., 2018; Mcsherry and Ritchie, 2013). However, between soil N and SOC was somehow discrete (r = 0.297; these effects were strongly context-specific, depending on p value = 0.001; R 2 = 0.088) in comparison to other stud- other factors like climate and soil-type vegetation (Abdalla ies (i.e. Yan et al., 2020). However, the novelty revealed by et al., 2018; Eze et al., 2018; Mcsherry and Ritchie, 2013). our model is that soil N could modulate the effects of certain Moreover, light and medium grazing intensities can increase SOC drivers, including livestock type and herbage NDF. SOC inputs by dung deposition and promoting aboveground Cattle-grazed grasslands stored more SOC than mixed and and root biomass production (Franzluebbers et al., 2000; sheep-grazed grasslands, but only under low soil N condi- Zeng et al., 2015). Considering that in our natural grasslands tions (Fig. 4b). Chang et al. (2018) found that in a N poor all grazing intensities are relatively low (see methods), our semi-arid grassland, sheep decreased SOC content in com- medium and high stock rates may increase soil carbon inputs parison to cattle due to vegetation changes caused by their Biogeosciences, 17, 6033–6050, 2020 https://doi.org/10.5194/bg-17-6033-2020 A. Rodríguez et al.: Soil organic carbon drivers 6045 feeding preference for highly palatable forbs (Sebastià et al., Our results concerning interactions between grazer type 2008a), thus promoting less palatable grasses which sup- and herbage quality provide some evidence of grazing ef- ported less root biomass. In overall, under low soil N con- fects not only through alterations of plant communities that ditions, palatable plants are expected to contribute to SOC were reported by previous studies in the region (Canals and inputs through the stimulation of C allocation in forb roots Sebastià, 2000; Sebastià et al., 2008a), but also through in- (Ågren and Franklin, 2003; Warembourg et al., 2003) and teractions with them. Although grazing effects were not the the increase in the overall plant productivity due to legume most important factors affecting SOC stocks, this is by far atmospheric N fixation (Van Der Heijden et al., 2008). the easiest component to manipulate in order to increase or However, these processes could decline under high soil maintain SOC in soils and face climate change (Komac et al., N contents. For instance, legume atmospheric N fixation 2014). Considering our results, we suggest conducting more could be reduced since it requires additional energy in com- experiments to investigate grazer-type effects on SOC under parison to nitrogen acquisition from the soil (Ibañez et al., different soil nutrient conditions and within plant commu- 2020; Minchin and Witty, 2005). Additionally, sheep selec- nities with contrasting herbage quality parameters. Grazing tive feeding habits could shift plant leaf traits in the com- management also has other advantages, such as preventing munity towards nutrient-conservative leaf traits, which com- the accumulation of aboveground C and reducing the risk of monly induce fungal-based soil food webs with slow nutrient forest fires (Nunes and Lourenço, 2017). cycling and high SOC storage due to low decomposition rates One key point of our results reinforces the idea that grazer (Orwin et al., 2010). type might be at least as important as grazing intensity in Additionally, grasslands with mixed grazed regimes stored regulating grassland ecosystem dynamics (Tóth et al., 2018) even more SOC than sheep-grazed grasslands under high soil and highlights the need for a more thorough research effort in N conditions (Fig. 4b, Table 3). This result emphasizes that disentangling not only grazing intensity, but also grazer-type mixed livestock assemblages deserve particular attention, be- effects on grassland soil organic carbon and nutrient cycling, cause mixed grazing can affect plant composition distinctly under different environmental circumstances. The combined from single grazing species regimes and alter traveling and model provided some evidence that grazing may affect SOC trampling behaviour of grazing animals (Aldezabal et al., not only through alterations of plant communities (Canals 2019; Chang et al., 2018; Liu et al., 2015). and Sebastià, 2000; Sebastià et al., 2008a), but also through NDF was negatively related to SOC at high soil N values interactions with them. Although grazing effects were not the (Fig. 4c). NDF proportion represents the amount of structural most important factors affecting SOC stocks, this is by far compounds on litter and hence is inversely related to non- the easiest component to manipulate in order to increase or structural compound content (Goering and Van Soest, 1970). maintain SOC in soils and face climate change (Komac et The latter are the main source of organic matter formation at al., 2014). Despite the need for precise knowledge on the ef- the early stages of decomposition, and they are incorporated fects of different land uses on ecosystems for climate change into microbial biomass in a highly efficient way (Cotrufo mitigation (Lo et al., 2015), studies addressing grazer-type et al., 2013). However, if microbial necromass was recycled effects on SOC are scarce (i.e. Zhou et al., 2017; Chang et by microbes before its incorporation into mineral-associated al., 2018). Considering our results, we suggest conducting organic matter (Córdova et al., 2018), it could be respired more experiments which investigate grazer-type effects on and mineralized instead of stored. Thus, our model suggests SOC under different soil nutrient conditions and within plant that incorporation of labile and fast metabolized non-organic communities with contrasting herbage quality parameters. 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 re- 5 Conclusion lated to SOC in sheep and mixed grazed grasslands (Fig. 4d). The ADL/NH ratio is a commonly used indicator for the re- The models presented here show a series of novel broad-scale sistance of litter to degradation, particularly at later stages of and local patterns concerning SOC stocks and their geo- decomposition (Taylor et al., 1989). Hence, the increase in physical, biochemical and grazing management drivers. Fac- SOC stocks with ADL/NH should be related to the physical tors driving SOC stocks often interacted in complex ways, pathway of soil organic matter incorporation, forming coarse within and between spatio-temporal scales. Temperature sea- particulate organic matter (Cotrufo et al., 2015). Moreover, sonality (TSIS) was the most critical geophysical factor, af- our model suggests that this pathway would be inhibited un- fecting SOC through interactions with topographical drivers der cattle grazing, presumably because of their higher diges- and grazing intensity. This illustrates how not only climate tive efficiency and thus less recalcitrant faeces (Wang et al., mean annual conditions should be considered when mod- 2018) and their less selective diet compared to sheep, as the elling SOC drivers, but also seasonal patterns. Concerning latter would avoid plants with a high lignin content, promot- biochemical factors, we found that the expected positive ing recalcitrant litter (Rosenthal et al., 2012; Sebastià et al., effect of soil N was modulated by livestock species and 2008a). herbage NDF, and herbage recalcitrance effects on SOC var- https://doi.org/10.5194/bg-17-6033-2020 Biogeosciences, 17, 6033–6050, 2020 6046 A. Rodríguez et al.: Soil organic carbon drivers ied depending on grazer type. Overall, we found a number of the Spanish Institute of Agronomical Research INIA (CARBO- interactions highlighting the need to expand knowledge on CLUS: SUM2006-00029-C02-0). Leticia San Emeterio was funded grassland SOC drivers under different conditions, especially through a Talent Recruitment grant from Obra Social La Caixa grazing. The latter is the most easily tractable factor affect- – Fundación CAN. The ARAID Foundation provided support to ing SOC. In conclusion, we provided valuable information Juan José Jiménez. This work was funded by the Spanish Science Foundation FECYT-MINECO (projects BIOGEI: GL2013-49142- for further studies dealing with SOC predictions at multi- C2-1-R and IMAGINE: CGL2017-85490-R) and the University of ple broad scales and laid out the basis for generating new Lleida (PhD Fellowship to Antonio Rodríguez). testable hypotheses for future studies, which may be useful for designing improved policies to palliate climate change. Review statement. This paper was edited by Frank Hagedorn and reviewed by two anonymous referees. Data availability. Data are not public as the PASTUS Database is currently being used for other research projects. Please contact one of us by e-mail for queries concerning the data used in this study. References Supplement. The supplement related to this article is available on- Abdalla, M., Hastings, A., Chadwick, D. R., Jones, D. L., Evans, line at: https://doi.org/10.5194/bg-17-6033-2020-supplement. C. D., Jones, M. B., Rees, R. M., and Smith, P.: Criti- cal review of the impacts of grazing intensity on soil or- ganic carbon storage and other soil quality indicators in exten- Author contributions. AR designed the statistical procedure, car- sively managed grasslands, Agr. Ecosyst. Environ., 253, 62–81, ried out the statistical analyses and wrote the original draft. RMC https://doi.org/10.1016/j.agee.2017.10.023, 2018. was responsible for field monitoring, lab analyses and acquisition Aerts, R. and Chapin, F. S.: The Mineral Nutrition of Wild Plants of information for the database implementation in the western Pyre- Revisited: A Re-evaluation of Processes and Patterns, Adv. Ecol. nees (Navarra). She also reviewed the draft. EA designed the NIRS Res., 30, 1–67, https://doi.org/10.1016/S0065-2504(08)60016-1, study and reviewed the draft. HD sampled and processed some of 1999. the data in the PASTUS Database and reviewed the draft. JGP pro- Ågren, G. I. and Franklin, O.: Root : shoot ratios, optimiza- cessed some of the data in the PASTUS Database and reviewed the tion and nitrogen productivity, Ann. Bot., 92, 795–800, draft. JP carried out the chemical analyses of herbage samples for https://doi.org/10.1093/aob/mcg203, 2003. NIR calibration and validation equations and reviewed the draft. Aldezabal, A., Pérez-López, U., Laskurain, N. A., and Odri- LSE designed the methodology and data collection and performed ozola, I.: El abandono del pastoreo afecta negativamente soil and vegetation sampling. She also reviewed the draft. ÀR sam- a la calidad del pasto en pastizales atlánticos ibéricos?; pled and processed some of the data in the PASTUS Database and Grazing abandonment negatively affects forage quality in reviewed the draft. JJJ collaborated in the fieldwork and reviewed iberian atlantic grasslands, Rev. Ecol. Montaña, 174, 42, the draft. MTS contributed to the conception, design and develop- https://doi.org/10.3989/pirineos.2019.174002, 2019. ment of the PASTUS Database. In addition, she ensured funding Bardgett, R. D. and Wardle, D. A.: Herbibore-mediated Linkages and coordinated the projects whose data are included in PASTUS. Between Abobeground and Belowground Communities, Ecol- Finally, she contributed to initial modelling, supervised the devel- ogy, 84, 2258–2268, https://doi.org/10.1890/02-0274, 2003. opment of the paper, and read and reviewed the drafts. Barton, K.: MuMIn: Multi-model inference, R package version 1.9.13, Version, 1.18, available at: https://cran.r-project.org/web/ packages/MuMIn/ (last access: 20 May 2019), 2015. Berninger, F., Susiluoto, S., Gianelle, D., Bahn, M., Wohlfahrt, G., Competing interests. The authors declare that they have no conflict Sutton, M., Garcia-Pausas, J., Gimeno, C., Sanz, M. J., Dore, S., of interest. Rogiers, N., Furger, M., Eugster, W., Balzarolo, M., Sebastià, M. T., Tenhunen, J., Staszewski, T., and Cernusca, A.: Management and site effects on carbon balances of european mountain mead- Acknowledgements. We would like to express our thanks to the ows and rangelands, Boreal Environ. Res., 20, 748–760, 2015. many people who collaborated in fieldwork, sample processing and Bing, H., Wu, Y., Zhou, J., Sun, H., Luo, J., Wang, J., and Yu, D.: data analysis over time. Stoichiometric variation of carbon, nitrogen, and phosphorus in soils and its implication for nutrient limitation in alpine ecosys- tem of Eastern Tibetan Plateau, J. Soils Sediments, 16, 405–416, Financial support. Research in this paper is based on the PAS- https://doi.org/10.1007/s11368-015-1200-9, 2016. TUS Database, which was compiled from different funding sources Bouyoucos, G. J.: Directions for making mechanical analysis of over time, the most relevant being the EU Interreg III-A Pro- soils by the hydrometer method, Soil Sci., 42, 225–2230, 1936. gramme (I3A-4-147-E) and the POCTEFA Programme/Interreg Burnham, K. P. and Anderson, D. R.: Model Selection and Mul- IV-A (FLUXPYR, EFA 34/08); the Spanish Science Foundation timodel Inference: A Practical Information-Theoretic Approach, FECYT-MICINN (CARBOPAS: REN2002-04300-C02-01; CAR- 2nd Edn., Springer, New York, 2002. BOAGROPAS: CGL2006-13555-C03-03 and CAPAS: CGL2010- Calcagno, V.: Package glmulti 1.0.7, Model selection and mul- 22378-C03-01); and the Foundation Catalunya-La Pedrera and timodel inference made easy, Community Ecol. Packag., 1– Biogeosciences, 17, 6033–6050, 2020 https://doi.org/10.5194/bg-17-6033-2020 A. Rodríguez et al.: Soil organic carbon drivers 6047 20, available at: https://cran.rproject.org/web/packages/glmulti/ bon sequestration potential, Sci. Total Environ., 565, 539–546, index.html (last access: 2 May 2019), 2015. https://doi.org/10.1016/j.scitotenv.2016.05.042, 2016. Cambardella, C. A. and Elliott, E. T.: Carbon and Nitrogen de Vries, F. T., Manning, P., Tallowin, J. R. B., Mortimer, S. Dynamics of Soil Organic Matter Fractions from Culti- R., Pilgrim, E. S., Harrison, K. A., Hobbs, P. J., Quirk, H., vated Grassland Soils, Soil Sci. Soc. Am. J., 58, 123–130, Shipley, B., Cornelissen, J. H. C., Kattge, J., and Bardgett, R. https://doi.org/10.2136/sssaj1994.03615995005800010017x, D.: Abiotic drivers and plant traits explain landscape-scale pat- 1994. terns in soil microbial communities, Ecol. Lett., 15, 1230–1239, Canals, R. M. and Sebastià, M. T.: Analyzing mechanisms regulat- https://doi.org/10.1111/j.1461-0248.2012.01844.x, 2012. ing diversity in rangelands through comparative studies: A case Díaz, S., Lavorel, S., De Bello, F., Quétier, F., Grigulis, K., and in the southwestern Pyrennees, Biodivers. Conserv., 9, 965–984, Robson, T. M.: Incorporating plant functional diversity effects https://doi.org/10.1023/A:1008967903169, 2000. in ecosystem service assessments, P. Natl. Acad. Sci. USA, 104, Carreras, J. and Diego, F. C.: Cartografia dels hàbitats de Catalunya 20684–20689, https://doi.org/10.1073/pnas.0704716104, 2007. 1 : 50 000, Generalitat de Catalunya, Barcelona, 2006. Duarte-guardia, S., Peri, P. L., Amelung, W., Sheil, D., Laffan, S. Chang, Q., Wang, L., Ding, S., Xu, T., Li, Z., Song, X., Zhao, X., W., Borchard, N., Bird, M. I., and Peri, P. L.: Better estimates of Wang, D., and Pan, D.: Grazer effects on soil carbon storage vary soil carbon from geographical data?: a revised global approach, by herbivore assemblage in a semi-arid grassland, J. Appl. Ecol., Mitig. Adapt. Strat. Gl., 24, 355–372, 2019. 55, 2517–2526, https://doi.org/10.1111/1365-2664.13166, 2018. Ellert, B. H. and Bettany, J. R.: Calculation of organic matter and Chang, R., Wang, G., Fei, R., Yang, Y., Luo, J., and Fan, J.: Alti- nutrients stored in soils under contrasting management regimes, tudinal Change in Distribution of Soil Carbon and Nitrogen in Can. J. Soil Sci., 75, 529–538, https://doi.org/10.4141/cjss95- Tibetan Montane Forests, Soil Sci. Soc. Am. J., 79, 1455–1469, 075, 1995. https://doi.org/10.2136/sssaj2015.02.0055, 2015. Elith, J., Leathwick, J. R., and Hastie, T.: A working guide Chapin, F. S., Bloom, A. J., Field, C. B., and Waring, R. H.: Plant to boosted regression trees, J. Anim. Ecol., 77, 802–813, Responses to Multiple Environmental Factors, Bioscience, 37, https://doi.org/10.1111/j.1365-2656.2008.01390.x, 2008. 49–57, https://doi.org/10.2307/1310177, 1987. Eze, S., Palmer, S. M., and Chapman, P. J.: Soil organic carbon Cleveland, C. C. and Liptzin, D.: C : N : P stoichiometry in soil: is stock in grasslands: Effects of inorganic fertilizers, liming and there a “Redfield ratio” for the microbial biomass?, Biogeochem- grazing in different climate settings, J. Environ. Manage., 223, istry, 85, 235–252, https://doi.org/10.1007/s10533-007-9132-0, 74–84, https://doi.org/10.1016/j.jenvman.2018.06.013, 2018. 2007. Fernández-Alonso, M. J., Díaz-Pinés, E., Ortiz, C., and Rubio, Córdova, S. C., Olk, D. C., Dietzel, R. N., Mueller, K. E., Archon- A.: Disentangling the effects of tree species and microclimate touilis, S. V., and Castellano, M. J.: Plant litter quality affects on heterotrophic and autotrophic soil respiration in a Mediter- the accumulation rate, composition, and stability of mineral- ranean ecotone forest, Forest Ecol. Manage., 430, 533–544, associated soil organic matter, Soil Biol. Biochem., 125, 115– https://doi.org/10.1016/j.foreco.2018.08.046, 2018. 124, https://doi.org/10.1016/J.SOILBIO.2018.07.010, 2018. Fernández-Martínez, M., Vicca, S., Janssens, I. A., Sardans, J., Cotrufo, M. F., Wallenstein, M. D., Boot, C. M., Denef, K., Luyssaert, S., Campioli, M., Chapin, F. S., Ciais, P., Malhi, Y., and Paul, E.: The Microbial Efficiency-Matrix Stabilization Obersteiner, M., Papale, D., Piao, S. L., Reichstein, M., Rodà, (MEMS) framework integrates plant litter decomposition with F., and Peñuelas, J.: Nutrient availability as the key regulator soil organic matter stabilization: do labile plant inputs form of global forest carbon balance, Nat. Clim. Chang., 4, 471–476, stable soil organic matter?, Glob. Change Biol., 19, 988–995, https://doi.org/10.1038/nclimate2177, 2014. https://doi.org/10.1111/gcb.12113, 2013. Fick, S. E. and Hijmans, R. J.: WorldClim 2: new 1-km spatial res- Cotrufo, M. F. F., Soong, J. L. J. L., Horton, A. J. A. J., Campbell, olution climate surfaces for global land areas, Int. J. Climatol., 3, E. E., Haddix, M. L. M. L. L., Wall, D. H. D. H., and Parton, 4302–4315, https://doi.org/10.1002/joc.5086, 2017. W. J. W. J.: Formation of soil organic matter via biochemical and Fox, J., Weisberg, S., Price, B., Adler, D., Bates, D., physical pathways of litter mass loss, Nat. Geosci., 8, 776–779, and Baud-Bovy, G.: Package “car,” R Doc., 146, https://doi.org/10.1038/ngeo2520, 2015. https://doi.org/10.1093/nar/gkq363, 2018. David, D. J.: The determination of exchangeable sodium, Franzluebbers, A. J., Stuedemann, J. A., Schomberg, H. H., and potassium, calcium and magnesium in soils by atomic- Wilkinson, S. R.: Soil organic C and N pools under long- absorption spectrophotometry, Analyst, 85, 495–503, term pasture management in the Southern Piedmont USA, Soil https://doi.org/10.1039/AN9608500495, 1960. Biol. Biochem., 32, 469–478, https://doi.org/10.1016/S0038- de Lamo, X. and Sebastià, M. T.: Role of biogeographical, ecolog- 0717(99)00176-5, 2000. ical and management factors on plant diversity in grasslands, in: Garcia-Pausas, J., Casals, P., Camarero, L., Huguet, C., Sebastià, Sociedad Española para el Estudio de los Pastos (SEEP), Vol. 11, M. T., Thompson, R., and Romanyà, J.: Soil organic car- edited by: Lloveras, J., González-Rodríguez, A., Vázquez-Yañez, bon storage in mountain grasslands of the Pyrenees: Effects O., Piñeiro, J., Santamaría, O., Olea, L., and Poblaciones, M. of climate and topography, Biogeochemistry, 82, 279–289, J., 832–834, Sociedad Española para el Estudio de los Pastos https://doi.org/10.1007/s10533-007-9071-9, 2007. (SEEP), Madrid, Spain, 2006. Garcia-Pausas, J., Romanyà, J., Montané, F., Rios, A. I., Taull, M., Deng, X., Zhan, Y., Wang, F., Ma, W., Ren, Z., Chen, X., Qin, Rovira, P., and Casals, P.: Are Soil Carbon Stocks in Moun- F., Long, W., Zhu, Z., and Lv, X.: Soil organic carbon of an tain Grasslands Compromised by Land-Use Changes?, in: High intensively reclaimed region in China: Current status and car- Mountain Conservation in a Changing World, edited by: Cata- https://doi.org/10.5194/bg-17-6033-2020 Biogeosciences, 17, 6033–6050, 2020 6048 A. Rodríguez et al.: Soil organic carbon drivers lan, J., Ninot, J. M., and Aniz, M. M., 207–230, Springer Inter- Appl., 10, 423–436, https://doi.org/10.1890/1051- national Publishing, Cham., 2017. 0761(2000)010[0423:TVDOSO]2.0.CO;2, 2000. Gobierno de Navarra: Mapa de usos del suelo de Navarra 1 : 25 000, Kennedy, M. J., Pevear, D. R., and Hill, R. J.: Mineral surface Gobierno de Navarra, Pamplona, Spain, 2003. control of organic carbon in black shale, Science, 295, 657–60, Goering, H. K. and Van Soest, P. J.: Forage fiber analyses (appa- https://doi.org/10.1126/science.1066611, 2002. ratus, reagents, prcedures, and some applications), USDA Agr Kikvidze, Z., Pugnaire, F. I., Brooker, R. W., Choler, P., Lortie, C. Handb, available at: http://agris.fao.org/agris-search/search.do? J., Michalet, R., and Callaway, R. M.: Linking patterns and pro- recordID=US201301156229 (last access: 18 April 2019), 1970. cesses in alpine plant communities: A global study, Ecology, 86, Gómez, D.: Aspectos ecológicos de los pastos, in: Pastos del Piri- 1395–1400, https://doi.org/10.1890/04-1926, 2005. neo, edited by: Fillat, F., García-González, R., Gómez, D., and Komac, B., Domènech, M., and Fanlo, R.: Effects of graz- Reiné, R., 61–73, Consejo Superior de Investigaciones Científi- ing on plant species diversity and pasture quality in sub- cas, CSIC, Madrid, Spain, 2008. alpine grasslands in the eastern Pyrenees (Andorra): Im- Grace, J. B.: Structural Equation Modeling and Natural Systems, plications for conservation, J. Nat. Conserv., 22, 247–255, Cambridge University Press, Cambridge, 2006. https://doi.org/10.1016/J.JNC.2014.01.005, 2014. Grace, J. B. and Bollen, K. A.: Interpreting the Results from Krull, E., Baldock, J., and Skjemstad, J.: Soil Texture Effects on Multiple Regression and Structural Equation Models, Bull. Decomposition and Soil Carbon Storage, in: Net Ecosystem Ecol. Soc. Am., 86, 283–295, https://doi.org/10.1890/0012- Exchange Workshop, 103–110, available at: http://citeseerx.ist. 9623(2005)86[283:itrfmr]2.0.co;2, 2005. psu.edu/viewdoc/summary?doi=10.1.1.121.3042 (last access: 13 Gray, J. M., Bishop, T. F. A., and Wilson, B. R.: Fac- May 2019), 2001. tors Controlling Soil Organic Carbon Stocks with Depth Lal, R.: Soil Carbon Sequestration Impacts on Global Climate in Eastern Australia, Soil Sci. Soc. Am. J., 79, 1741, Change and Food Secruity, Science, 304, 1623–1627, 2004a. https://doi.org/10.2136/sssaj2015.06.0224, 2015. Lal, R.: Soil carbon sequestration to miti- Greenwell, B., Boehmke, B., and Cunningham, J.: Package “gbm”, gate climate change, Geoderma, 123, 1–22, available at: https://cran.r-project.org/web/packages/gbm/index. https://doi.org/10.1016/J.GEODERMA.2004.01.032, 2004b. html, last access: 15 February 2019. Legendre, P. and Legendre, L.: Numerical Ecology, 2nd Edn., Else- Grömping, U.: R package relaimpo: relative importance vier, Amsterdam, 1998. for linear regression, J. Stat. Softw., 17, 139–147, Leifeld, J., Meyer, S., Budge, K., Sebastia, M. T., Zimmer- https://doi.org/10.1016/j.foreco.2006.08.245, 2006. mann, M., and Fuhrer, J.: Turnover of grassland roots in Grueber, C. E., Nakagawa, S., Laws, R. J., and Jamieson, I. G.: Mul- mountain ecosystems revealed by their radiocarbon signature: timodel inference in ecology and evolution: Challenges and solu- Role of temperature and management, PLoS One, 10, 1–13, tions, J. Evol. Biol., 24, 699–711, https://doi.org/10.1111/j.1420- https://doi.org/10.1371/journal.pone.0119184, 2015. 9101.2010.02210.x, 2011. Lenth, R., Singmann, H., Love, J., Buerkne, P., and Herve, M.: Haden, A. C., Yang, W. H., and DeLucia, E. H.: Soils’ dirty Package “emmeans”, available at: https://cran.r-project.org/web/ little secret: Depth-based comparisons can be inadequate for packages/emmeans/, last access: 18 March 2019. quantifying changes in soil organic carbon and other min- Liu, J., Feng, C., Wang, D., Wang, L., Wilsey, B. J., and Zhong, eral soil properties, Glob. Change Biol., 26, 3759–3770, Z.: Impacts of grazing by different large herbivores in grassland https://doi.org/10.1111/gcb.15124, 2020. depend on plant species diversity, J. Appl. Ecol., 52, 1053–1062, Hijmans, R. J., Phillips, S., Leathwick, J., and Maintainer, J. E.: https://doi.org/10.1111/1365-2664.12456, 2015. Package “dismo” Type Package, Species Distribution Model- Lo, Y. H., Blanco, J. A., Canals, R. M., González de Andrés, E., San ing, available at: https://cran.r-project.org/web/packages/dismo/ Emeterio, L., Imbert, J. B., and Castillo, F. J.: Land use change dismo.pdf (lsat access: 15 February 2019), 2017. effects on carbon and nitrogen stocks in the Pyrenees during the Hobley, E., Wilson, B., Wilkie, A., Gray, J., and Koen, T.: last 150 years: A modeling approach, Ecol. Modell., 312, 322– Drivers of soil organic carbon storage and vertical dis- 334, https://doi.org/10.1016/j.ecolmodel.2015.06.005, 2015. tribution in Eastern Australia, Plant Soil, 390, 111–127, López-Moreno, J. I., Revuelto, J., Gilaberte, M., Morán-Tejeda, E., https://doi.org/10.1007/s11104-015-2380-1, 2015. Pons, M., Jover, E., Esteban, P., García, C., and Pomeroy, J. W.: Hook, P. B. and Burke, I. C.: Biogeochemistry in a shortgrass The effect of slope aspect on the response of snowpack to climate landscape: control by topography, soil texture, and microcli- warming in the Pyrenees, Theor. Appl. Climatol., 117, 207–219, mate, Ecology, 81, 2686–2703, https://doi.org/10.1890/0012- https://doi.org/10.1007/s00704-013-0991-0, 2013. 9658(2000)081[2686:BIASLC]2.0.CO;2, 2000. Lozano-García, B., Parras-Alcántara, L., and Brevik, E. C.: Im- Ibañez, M., Altimir, N., Ribas, A., Eugster, W., and Sebastià, M. pact of topographic aspect and vegetation (native and refor- T.: Phenology and plant functional type dominance drive CO2 ested areas) on soil organic carbon and nitrogen budgets in exchange in seminatural grasslands in the Pyrenees, J. Agr. Sci., Mediterranean natural areas, Sci. Total Environ., 544, 963–970, 158, 3–14, https://doi.org/10.1017/S0021859620000179, 2020. https://doi.org/10.1016/J.SCITOTENV.2015.12.022, 2016. Jenny, H.: Factors of soil formation: A System of Quantitative Lu, X., Yan, Y., Sun, J., Zhang, X., Chen, Y., Wang, X., and Cheng, Pedology, Dover Publications, New York, USA, 1941. G.: Carbon, nitrogen, and phosphorus storage in alpine grassland Jobbágy, E. G. and Jackson, R. B.: The Vertical Dis- ecosystems of Tibet: effects of grazing exclusion, Ecol. Evol., 5, tribution of Soil Organic Carbon and Its Relation 4492–504, https://doi.org/10.1002/ece3.1732, 2015. to Climate and Vegetation, Source Ecol. Appl. Ecol. Lu, X., Kelsey, K. C., Yan, Y., Sun, J., Wang, X., Cheng, G., and Neff, J. C.: Effects of grazing on ecosystem structure and func- Biogeosciences, 17, 6033–6050, 2020 https://doi.org/10.5194/bg-17-6033-2020 A. Rodríguez et al.: Soil organic carbon drivers 6049 tion of alpine grasslands in Qinghai-Tibetan Plateau: a synthesis, Rosenthal, G., Schrautzer, J., and Eichberg, C.: Low-intensity graz- Ecosphere, 8, e01656, https://doi.org/10.1002/ecs2.1656, 2017. ing with domestic herbivores: A tool for maintaining and restor- Manning, P., de Vries, F. T., Tallowin, J. R. B., Smith, R., Mor- ing plant diversity in temperate Europe, Tuexenia, 32, 167–205, timer, S. R., Pilgrim, E. S., Harrison, K. A., Wright, D. G., 2012. Quirk, H., Benson, J., Shipley, B., Cornelissen, J. H. C., Kattge, Sanderman, J., Amundson, R. G., and Baldocchi, D. D.: Application J., Bönisch, G., Wirth, C., and Bardgett, R. D.: Simple mea- of eddy covariance measurements to the temperature dependence sures of climate, soil properties and plant traits predict national- of soil organic matter mean residence time, Global Biogeochem. scale grassland soil carbon stocks, J. Appl. Ecol., 52, 1188–1196, Cy., 17, 1061, https://doi.org/10.1029/2001GB001833, 2003. https://doi.org/10.1111/1365-2664.12478, 2015. Schlesinger, W. H.: Carbon Balance in Terres- Mcsherry, M. E. and Ritchie, M. E.: Effects of grazing on grassland trial Detritus, Annu. Rev. Ecol. Syst., 8, 51–81, soil carbon: A global review, Glob. Change Biol., 19, 1347–1357, https://doi.org/10.1146/annurev.es.08.110177.000411, 1977. https://doi.org/10.1111/gcb.12144, 2013. Schöning, I., Grüneberg, E., Sierra, C. A., Hessenmöller, D., Minchin, F. R. and Witty, J. F.: Respiratory/carbon costs of symbi- Schrumpf, M., Weisser, W. W., and Schulze, E. D.: Causes of otic nitrogen fixation in legumes, in: Plant respiration, 195–205, variation in mineral soil C content and turnover in differently Springer, Dordrecht, the Netherlands, 2005. managed beech dominated forests, Plant Soil, 370, 625–639, Nunes, A. N. and Lourenço, L.: Increased vulnerability to wild- https://doi.org/10.1007/s11104-013-1654-8, 2013. fires and post fire hydro-geomorphic processes in Portuguese Sebastià, M.-T.: Role of topography and soils in grassland structur- mountain regions: What has changed?, Open Agric., 2, 70–82, ing at the landscape and community scales, Basic Appl. Ecol., 5, https://doi.org/10.1515/opag-2017-0008, 2017. 331–346, https://doi.org/10.1016/j.baae.2003.10.001, 2004. Olsen, S. R.: Estimation of available phosphorus in soils by extrac- Sebastià, M.-T., Kirwan, L., and Connolly, J.: Strong shifts tion with sodium bicarbonate, US Dept. of Agriculture, Wash- in plant diversity and vegetation composition in grassland ington D.C., USA, 1954. shortly after climatic change, J. Veg. Sci., 19, 299–306, Orwin, K. H., Buckland, S. M., Johnson, D., Turner, B. L., https://doi.org/10.3170/2008-8-18356, 2008a. Smart, S., Oakley, S., and Bardgett, R. D.: Linkages of plant Sebastià, M.-T., de Bello, F., Puig, L., and Taull, M.: Grazing as a traits to soil properties and the functioning of temperate grass- factor structuring grasslands in the Pyrenees, Appl. Veg. Sci., 11, land, J. Ecol., 98, 1074–1083, https://doi.org/10.1111/j.1365- 215–222, https://doi.org/10.3170/2008-7-18358, 2008b. 2745.2010.01679.x, 2010. Shipley, B. and Parent, M.: Germination Responses of 64 Wetland Ottoy, S., Van Meerbeek, K., Sindayihebura, A., Hermy, M., and Species in Relation to Seed Size, Minimum Time to Reproduc- Van Orshoven, J.: Assessing top- and subsoil organic carbon tion and Seedling Relative Growth Rate, Funct. Ecol., 5, 111, stocks of Low-Input High-Diversity systems using soil and https://doi.org/10.2307/2389561, 1991. vegetation characteristics, Sci. Total Environ., 589, 153–164, Simard, R. R.: Ammonium acetate-extractable elements, in Soil https://doi.org/10.1016/J.SCITOTENV.2017.02.116, 2017. sampling and methods of analysis, Vol. 1, 39–42, Lewis Pub- Penin, D.: Cartographie des habitats naturels (typologie Corine lisher FL, USA, 1993. Biotope) du massif Madres-Coronat, 1st Edn., AGRNN/CPIE du Solly, E. F., Schöning, I., Boch, S., Kandeler, E., Marhan, S., Conflent, AGRNN/CPIE du Conflent Paris, France, 1997. Michalzik, B., Müller, J., Zscheischler, J., Trumbore, S. E., and Peri, P. L., Rosas, Y. M., Ladd, B., Toledo, S., Lasagno, R. G., Schrumpf, M.: Factors controlling decomposition rates of fine and Pastur, G. M.: Modelling soil carbon content in South root litter in temperate forests and grasslands, Plant Soil, 382, Patagonia and evaluating changes according to climate, veg- 203–218, https://doi.org/10.1007/s11104-014-2151-4, 2014. etation, desertification and grazing, Sustainability, 10, 438, Stockmann, U., Adams, M. A., Crawford, J. W., Field, D. J., https://doi.org/10.3390/su10020438, 2018. Henakaarchchi, N., Jenkins, M., Minasny, B., McBratney, A. Puissant, J., Jassey, V. E. J., Mills, R. T. E., Robroek, B. J. M., B., Courcelles, V. de R. de, Singh, K., Wheeler, I., Abbott, Gavazov, K., De Danieli, S., Spiegelberger, T., Griffiths, R., L., Angers, D. A., Baldock, J., Bird, M., Brookes, P. C., Buttler, A., Brun, J.-J. J., and Cécillon, L.: Seasonality alters Chenu, C., Jastrow, J. D., Lal, R., Lehmann, J., O’Donnell, drivers of soil enzyme activity in subalpine grassland soil un- A. G., Parton, W. J., Whitehead, D., and Zimmermann, M.: dergoing climate change, Soil Biol. Biochem., 124, 266–274, The knowns, known unknowns and unknowns of sequestra- https://doi.org/10.1016/j.soilbio.2018.06.023, 2018. tion of soil organic carbon, Agr. Ecosyst. Environ., 164, 80–99, R Core Team: R: A Language and Environment for Statistical https://doi.org/10.1016/j.agee.2012.10.001, 2013. Computing, available at: https://www.r-project.org/, last access: Symonds, M. R. E. and Moussalli, A.: A brief guide to model selec- 15 February 2016. tion, multimodel inference and model averaging in behavioural Rillig, M. C., Ryo, M., Lehmann, A., Aguilar-Trigueros, C. A., ecology using Akaike’s information criterion, Behav. Ecol. So- Buchert, S., Wulf, A., Iwasaki, A., Roy, J., and Yang, G.: ciobiol., 65, 13–21, https://doi.org/10.1007/s00265-010-1037-6, The role of multiple global change factors in driving soil 2011. functions and microbial biodiversity, Science, 366, 886–890, Taylor, B. R., Parkinson, D., and Parsons, W. F. J.: Ni- https://doi.org/10.1126/science.aay2832, 2019. trogen and Lignin Content as Predictors of Litter De- Rodríguez, A., de Lamo, X., Sebastiá, M.-T., and Sebastià, M.- cay Rates: A Microcosm Test, Ecology, 70, 97–104, T.: Interactions between global change components drive plant https://doi.org/10.2307/1938416, 1989. species richness patterns within communities in mountain grass- Tipping, E., Somerville, C. J., and Luster, J.: The C : N : P : S sto- lands independetly of topography, J. Veg. Sci., 29, 1029–1039, ichiometry of soil organic matter, Biogeochemistry, 130, 117– https://doi.org/10.1111/jvs.12683, 2018. 131, https://doi.org/10.1007/s10533-016-0247-z, 2016. https://doi.org/10.5194/bg-17-6033-2020 Biogeosciences, 17, 6033–6050, 2020 6050 A. Rodríguez et al.: Soil organic carbon drivers Tóth, E., Deák, B., Valkó, O., Kelemen, A., Miglécz, T., Tóth- data-assimilation analysis, Geoderma, 262, 235–242, mérész, B., and Török, P.: Livestock Type is More Crucial https://doi.org/10.1016/j.geoderma.2015.08.038, 2016a. Than Grazing Intensity: Traditional Cattle and Sheep Graz- Xu, X., Shi, Z., Li, D., Rey, A., Ruan, H., Craine, J. M., ing in Short-Grass Steppes, Land Degrad. Dev., 29, 231–239, Liang, J., Zhou, J., and Luo, Y.: Soil properties con- https://doi.org/10.1002/ldr.2514, 2018. trol decomposition of soil organic carbon: Results from Van Der Heijden, M. G. A., Bardgett, R. D., and Van Straalen, data-assimilation analysis, Geoderma, 262, 235–242, N. M.: The unseen majority: Soil microbes as drivers of plant https://doi.org/10.1016/J.GEODERMA.2015.08.038, 2016b. diversity and productivity in terrestrial ecosystems, Ecol. Lett., Yan, J., Wang, L., Hu, Y., Tsang, Y. F., Zhang, Y., Wu, J., Fu, X., and 11, 296–310, https://doi.org/10.1111/j.1461-0248.2007.01139.x, Sun, Y.: Plant litter composition selects different soil microbial 2008. structures and in turn drives different litter decomposition pattern Vitousek, P. M. and Howarth, R. W.: Nitrogen Limitation on Land and soil carbon sequestration capability, Geoderma, 319, 194– and in the Sea: How Can It Occur?, 1991. 203, https://doi.org/10.1016/J.GEODERMA.2018.01.009, 2018. Von Haden, A. C. and Dornbush, M. E.: Patterns of root Yan, L., Li, Y., Wang, L., Zhang, X., Wang, J., Wu, H., Yan, Z., decomposition in response to soil moisture best explain Zhang, K., and Kang, X.: Grazing significantly increases root high soil organic carbon heterogeneity within a mesic, shoot ratio but decreases soil organic carbon in Qinghai-Tibetan restored prairie, Agr. Ecosyst. Environ., 185, 188–196, Plateau grasslands: A hierarchical meta-analysis, Land Degrad. https://doi.org/10.1016/j.agee.2013.12.027, 2014. Dev., 31, 2369–2378, https://doi.org/10.1002/ldr.3606, 2020. Van Soest, P. J., Robertson, J. B., and Lewis, B. A.: Meth- Yan, Y., Tian, L., Du, Z., Chang, S. X., and Cai, Y.: Carbon, nitro- ods for Dietary Fiber, Neutral Detergent Fiber, and Non- gen and phosphorus stocks differ among vegetation patch types starch Polysaccharides in Relation to Animal Nutrition, J in a degraded alpine steppe, J. Soils Sediments, 19, 1809–1819, Dairy Sci., 74, 3583–3597, https://doi.org/10.3168/jds.S0022- https://doi.org/10.1007/s11368-018-2191-0, 2019. 0302(91)78551-2, 1991. Yang, Y., Chen, Y., Li, Z., and Chen, Y.: Land-use/cover conver- Wang, J., Wang, D., Li, C., Seastedt, T. R., Liang, C., Wang, L., sion affects soil organic-carbon stocks: A case study along the Sun, W., Liang, M., and Li, Y.: Feces nitrogen release induced main channel of the Tarim River, China, PLoS One, 13, 1–14, by different large herbivores in a dry grassland, Ecol. Appl., 28, https://doi.org/10.1371/journal.pone.0206903, 2018. 201–211, https://doi.org/10.1002/eap.1640, 2018. Yang, Y. H., Fang, J. Y., Pan, Y. D., and Ji, C. J.: Aboveground Wanyama, I., Pelster, D. E., Butterbach-Bahl, K., Verchot, L. biomass in Tibetan grasslands, J. Arid Environ., 73, 91–95, V, Martius, C., and Rufino, M. C.: Soil carbon dioxide and https://doi.org/10.1016/J.JARIDENV.2008.09.027, 2009. methane fluxes from forests and other land use types in an Yuan, Z.-Q., Fang, C., Zhang, R., Li, F.-M., Javaid, M. African tropical montane region, Biogeochemistry, 143, 171– M., and Janssens, I. A.: Topographic influences on soil 190, https://doi.org/10.1007/s10533-019-00555-8, 2019. properties and aboveground biomass in lucerne-rich vegeta- Warembourg, F. R., Roumet, C., and Lafont, F.: Differ- tion in a semi-arid environment, Geoderma, 344, 137–143, ences in rhizosphere carbon-partitioning among plant https://doi.org/10.1016/J.GEODERMA.2019.03.003, 2019. species of different families, Plant Soil, 256, 347–357, Zeng, C., Wu, J., and Zhang, X.: Effects of grazing on https://doi.org/10.1023/A:1026147622800, 2003. above-vs. below-ground biomass allocation of alpine grass- White, R., Murray, S., and Rohweder, M.: Pilot Analysis of Global lands on the northern tibetan plateau, PLoS One, 10, 1–15, Ecosystems: Grassland Ecosystems, 2000. https://doi.org/10.1371/journal.pone.0135173, 2015. Wiesmeier, M., Urbanski, L., Hobley, E., Lang, B., von Lützow, Zhang, X., Liu, M., Zhao, X., Li, Y., Zhao, W., Li, A., Chen, M., Marin-Spiotta, E., van Wesemael, B., Rabot, E., Ließ, M., S., Chen, S., Han, X., and Huang, J.: Topography and graz- Garcia-Franco, N., Wollschläger, U., Vogel, H. J., and Kögel- ing effects on storage of soil organic carbon and nitrogen Knabner, I.: Soil organic carbon storage as a key function of soils in the northern China grasslands, Ecol. Indic., 93, 45–53, – A review of drivers and indicators at various scales, Geoderma, https://doi.org/10.1016/J.ECOLIND.2018.04.068, 2018. 333, 149–162, https://doi.org/10.1016/j.geoderma.2018.07.026, Zhao, Y., Ding, Y., Hou, X., Li, F. Y., Han, W., and Yun, X.: Ef- 2019. fects of temperature and grazing on soil organic carbon storage in Williams, P. C. and Sobering, D. C.: How do we do it: a brief sum- grasslands along the Eurasian steppe eastern transect, PLoS One, mary of the methods we use in developing near infrared calibra- 12, 1–16, https://doi.org/10.1371/journal.pone.0186980, 2017. tions, in Near infrared spectroscopy: The future waves, edited Zhou, G., Zhou, X., He, Y., Shao, J., Hu, Z., Liu, R., Zhou, H., by: Williams, P. C., Davies, A. M. C., and Williams, P., 185–188, and Hosseinibai, S.: Grazing intensity significantly affects be- NIR Publications, Chichester, UK, 1996. lowground carbon and nitrogen cycling in grassland ecosys- Williams, T. E., Chang, C. M., Rosen, E. L., Garcia, G., Runner- tems: a meta-analysis, Glob. Change Biol., 23, 1167–1179, strom, E. L., Williams, B. L., Koo, B., Buonsanti, R., Milliron, https://doi.org/10.1111/gcb.13431, 2017. D. J., and Helms, B. A.: NIR-Selective electrochromic hetero- Zhu, M., Feng, Q., Qin, Y., Cao, J., Zhang, M., Liu, W., Deo, R. C., material frameworks: A platform to understand mesoscale trans- Zhang, C., Li, R., and Li, B.: The role of topography in shaping port phenomena in solid-state electrochemical devices, J. Mater. the spatial patterns of soil organic carbon, Catena, 176, 296–305, Chem. C, 2, 3328–3335, https://doi.org/10.1039/c3tc32247e, https://doi.org/10.1016/J.CATENA.2019.01.029, 2019. 2014. Xu, X., Shi, Z., Li, D., Rey, A., Ruan, H., Craine, J. M., Liang, J., Zhou, J., and Luo, Y.: Soil properties con- trol decomposition of soil organic carbon: Results from Biogeosciences, 17, 6033–6050, 2020 https://doi.org/10.5194/bg-17-6033-2020