Improved management increases soil mineral-protected organic carbon storage via plant-microbial-nutrient mediation in semi-arid grasslands

Soil organic carbon (SOC) storage in semi-arid grasslands is threatened by both climate change and land degradation, impacting food production and climate regulation. Improved management has been proposed to increase SOC stocks and overcome these challenges. However, the benefits of improved management practices in semi-arid regions are in question. Little is known about the effects of management on the functional components of SOC, particulate (POC) and mineral-associated organic carbon (MAOC), which are expected to respond differently, and about the pathways that mediate these responses, such as changes in vegetation and soil microbial communities.

This work analyses the effect of rotational grazing, legume sowing and grazing exclusion on topsoil (0–8 cm) SOC, POC and MAOC stocks in Mediterranean wooded grasslands compared to continuous conventional grazing. Changes in plant diversity and morpho-biochemical traits, soil fertility and microbial composition were also evaluated. A total of 188 plots were sampled in 9 farms across a wide environmental gradient.

More resource-acquisitive, nitrogen-rich and less lignified plant community, higher soil microbial biomass with lower Gram+ $/$ Gram− ratio, and higher soil fertility were associated with higher SOC storage, with similar impacts on POC and MAOC. Rotational grazing increased POC, MAOC and total SOC stocks by 15 %, 15 % and 14 % respectively, compared to continuous grazing. This effect was mediated by an increase in soil fertility in the rotationally grazed paddocks. On the other hand, grazing exclusion reduced POC stocks by 12 % compared to continuous grazing. This depletion was mainly due to a reduction in microbial biomass and litter quality of vegetation in non-grazed paddocks. Both POC and MAOC stocks were lower at the warmer sites.

We conclude that rotational grazing can enhance long-term SOC storage in semi-arid grasslands, thereby increasing their resilience and climate mitigation capacity, whereas abandoning grazing could lead to SOC losses.

1 Introduction

Grasslands cover 40.5 % of the world's ice-free land area and store one third of the terrestrial carbon (C) stocks, mainly (over 80 %) in the form of soil organic matter (SOM) (White et al., 2000). This C pool in grassland soils surpasses the global aboveground vegetation C pool (IPCC, 2023a), highlighting the importance of grasslands SOM in global climate regulation (Bai and Cotrufo, 2022). In addition, SOM is a cornerstone of grassland productivity and functioning (Cotrufo and Lavallee, 2022; Tiessen et al., 1994). However, C storage in grassland soils is threatened by increasing C outputs (i.e., soil respiration and erosion) and decreasing C inputs (i.e., primary productivity) due to anthropogenic impacts such as land degradation and climate change (Crowther et al., 2016; Gang et al., 2014; Godde et al., 2020; Lei et al., 2021). Improved grassland management offers a key opportunity to counter these threats by increasing soil organic C (SOC) stocks and supporting both climate change mitigation and adaptation (Bai and Cotrufo, 2022; Conant, 2012; Dondini et al., 2023; Stanley et al., 2024). This is particularly valuable in the case of semi-arid grasslands, which represent the majority of the global grassland area (White et al., 2000) but are also more vulnerable to climate change impacts than wetter grasslands (Smith et al., 2024).

Practices such as rotational grazing, legume sowing and grazing exclusion have been proposed to increase and conserve SOC stocks (Conant et al., 2017; Yu et al., 2021). Rotational grazing encompasses a variety of practices where small paddocks are grazed at high intensity for short periods of time, allowing for longer pasture rest than continuous grazing (Teague et al., 2013). This management favors vegetation recovery after defoliation and reduces grazing patchiness and livestock selectivity (Jacobo et al., 2006; Teague et al., 2013). The effects of rotational grazing on grassland productivity and animal performance are under debate (Briske et al., 2008; Teague et al., 2013; di Virgilio et al., 2019) but most studies have found positive effects on SOC stocks associated with this practice (Byrnes et al., 2018; Conant et al., 2017; Phukubye et al., 2022; Teague et al., 2011). Sowing legumes in natural pastures is practiced worldwide, having clear positive impacts on SOC stocks (Carranca et al., 2022; Conant et al., 2017; Moreno et al., 2021) and grassland productivity (Bartholomew and Williams, 2010; Carrascosa et al., 2024; Jaurena et al., 2016; Khatiwada et al., 2020; Rama et al., 2022). Grazing exclusion is widely advocated as a tool for ecosystem restoration (Cheng et al., 2016; Novelly and Watson, 2007) and is also a global trend driven by land abandonment, particularly in high-income countries (Li and Li, 2017). The effects of grazing exclusion on SOC in grasslands are mixed, showing both positive (Cheng et al., 2016; Yu et al., 2021) and negative outcomes (Wilson et al., 2018). Moreover, the net effect of management practices has been shown to depend on the environmental context (Maestre et al., 2022; McSherry and Ritchie, 2013; Niu et al., 2025), and global meta-analyses remain inconclusive. While some global studies have reported greater benefits on SOC stocks from grazing exclusion and rotational grazing in arid climates (Zhou et al., 2017), others found these management practices to be more beneficial in wetter climates (Byrnes et al., 2018; Hawkins, 2017; Zhou et al., 2019). To the best of our knowledge, no similar studies have addressed the interaction between environmental conditions and the effect of legume sowings on SOC stocks. Thus, some questions remain open and further research is needed to clarify the net effects of improved management practices on SOC storage under different environmental conditions and particularly in semi-arid grasslands.

To better understand the controls and vulnerability of C in soils, SOM can be conceptualized into particulate organic matter (POM) and mineral-associated organic matter (MAOM) fractions (Cotrufo and Lavallee, 2022; Lavallee et al., 2020). POM originates from fragmented structural plant inputs and, to a lesser extent, microbial recalcitrant compounds (Angst et al., 2021; Six et al., 2001). In contrast, MAOM forms through the sorption of microbial necromass and soluble plant inputs onto soil mineral surfaces (Angst et al., 2021; Cotrufo et al., 2022). The mineral bonds partially protect MAOM from decomposition (Baldock and Skjemstad, 2000), while POM is readily accessible to microbial degradation, although occlusion within soil aggregates can reduce its accessibility (Angst et al., 2017). Hence, POM accumulation is mainly controlled by environmental constraints on microbial activity, e.g. low temperatures and highly acidic soils (Hansen et al., 2024; Yu et al., 2022; Zhou et al., 2024). For the same reasons, POM is more vulnerable to climatic warming (Benbi et al., 2014; Georgiou et al., 2024; Rocci et al., 2021) and the mean residence time of C in POM (years to decades) is on average shorter than in MAOM (decades to centuries; Zhou et al., 2024). However, MAOM storage in soil is theoretically limited by the availability of free mineral surface area (i.e. clay and silt content; Six et al., 2002), and a saturation point can be observed, where no more MAOM accumulate despite increases in total SOM contents (Cotrufo et al., 2019; Georgiou et al., 2022). Management changes have been shown to influence SOM fractions in grassland soils (Khatri-Chhetri et al., 2024; Mosier et al., 2021; Oliveira Filho et al., 2019), and a conceptual framework for relating grazing management to SOM distribution has recently been proposed (Stanley et al., 2024). Yet, the underlying processes mediating these effects, such as alterations in vegetation or soil microbial communities (Laliberté and Tylianakis, 2012; Peco et al., 2017; Wilson et al., 2018), remain poorly understood.

Primary productivity is the point of entry of C into soil, and consequently the amount of plant inputs regulates SOM accrual (King et al., 2023; Zhou et al., 2024), but microbial processing largely determines the fate of that C (Cotrufo and Lavallee, 2022; Crowther et al., 2019). In this sense, the chemical composition of plant inputs, and the soil microbiota carbon use efficiency (CUE), i.e. the amount of C used for microbial growth and products relative to total C uptake, play a crucial role in the SOM formation process (Cotrufo and Lavallee, 2022; Tao et al., 2023). Recalcitrant plant inputs [i.e., high carbon-nitrogen ratio (C $/$ N) and lignin content] tend to promote short-term SOM accumulation, primarily as POM, due to their chemical resistance to decomposition (Cheng et al., 2023; Cotrufo and Lavallee, 2022). However, as outlined in the Microbial Efficiency-Matrix Stabilization (MEMS) framework (Cotrufo et al., 2013), recalcitrant inputs are less efficiently decomposed by microbes, leading to greater C losses in the long term, compared to labile (i.e., water-soluble, low C $/$ N and lignin content) plant inputs (Cotrufo and Lavallee, 2022; Ridgeway et al., 2022). Thus, labile plant inputs are expected to enhance MAOM formation and SOM stocks in the long term, due to their faster and more efficient decomposition (Cheng et al., 2023; Haddix et al., 2016). Elias et al. (2024) added complexity to these assumptions, showing that plant input characteristics may favor certain microbial groups over others, altering the overall CUE of the microbial community. For example, fungi, which are often assumed to have a higher CUE than bacteria (Kallenbach et al., 2016; Strickland and Rousk, 2010), are favored by the addition of recalcitrant inputs (Bai et al., 2024; Strickland and Rousk, 2010). Substrate preferences have been also identify for Gram-positive (Gram+) and Gram-negative (Gram−) bacteria (Fanin et al., 2019; Kramer and Gleixner, 2008), with consequences for SOC accrual (Klumpp et al., 2009). Importantly, much of the research on the influence of plant input characteristics and microbial communities on SOM formation dynamics has relied on incubation experiments (Cheng et al., 2023; Haddix et al., 2016; Ridgeway et al., 2022) and there is limited information on how these findings translate to natural field conditions. Other vegetation characteristics such as species richness have been shown to positively influence SOC stocks (Lange et al., 2015; Steinbeiss et al., 2008), but its effects on SOC fractions have been poorly evaluated, with inconclusive results in grasslands (Mortensen et al., 2025). In addition, the relationships between SOM stocks and fractions and plant functional traits have rarely been studied (Manning et al., 2015; Mortensen et al., 2025; Xu et al., 2021) despite the latter being widely used to predict ecosystem functioning and responses (Funk et al., 2017). Plant functional traits are highly correlated with processes such as litter decomposition (Cornwell et al., 2008; Fortunel et al., 2009; Kazakou et al., 2009) or root exudates production (Guyonnet et al., 2018) and may therefore be a promising tool to study the relationships between vegetation, soil microbiota and SOM formation dynamics (Faucon et al., 2017).

The aim of this work is to evaluate the impact of rotational grazing, legume sowing, and grazing exclusion on topsoil (0–8 cm) SOC stocks, and their distribution between POM and MAOM fractions, in semi-arid grasslands, compared to conventional continuous grazing. We also evaluate changes in vegetation characteristics (identity and diversity of biochemical and morphological traits), soil nutrients, and soil microbial communities as possible pathways through which management might indirectly affect SOC stocks and fractions. This study focuses on the Iberian dehesas Mediterranean woody grassland, the main example of semi-arid grasslands in Europe (Porqueddu et al., 2016). The most widespread livestock management in this ecosystem is continuous grazing, but in recent decades rotational grazing and legume sowing have gained importance (Frongia et al., 2023; Pulina et al., 2023). At the same time, the amount of ungrazed pastures has increased due to land abandonment (Palomo-Campesino et al., 2018). Iberian dehesas occupy 3.1 million ha, spanning a wide environmental gradient, and have been subject to extensive grazing for centuries (Moreno and Pulido, 2009), making them an ideal model system for assessing the effects of improved management on SOC stocks in semi-arid grasslands. In particular, we designed our study to answer the following questions: (1) what are the effects of the different management practices on bulk SOC and fractions stocks? (2) Are these effects mediated by changes in vegetation or soil microbial communities? (3) Is SOC storage in these grasslands, and its enhancement, modulated or limited by environmental factors such as climate or soil properties? And (4) are these mechanisms and controls the same or different for C in POM or MAOM? Understanding the potential and limitations of improved management on SOC storage in POM and MAOM in semi-arid grasslands can guide policymakers in enhancing the climate change adaptation and mitigation capacity of these ecosystems, while supporting productivity and soil fertility.

2 Methods

2.1 Study area and experimental design

The study was carried out at nine commercial dehesa farms located along a north–south gradient in the western part of the Iberian Peninsula (Fig. 1a). The region has a continental Mediterranean climate, but on a local scale, in relative terms, farms can be grouped into three main climatic regions (Fig. 1b and c) according to the average climate for the period 1980–2018 (García Bravo et al., 2023). A cold-dry region – 12.9 °C mean annual temperature (MAT) and 445 mm mean annual precipitation (MAP) – in the north; a warm-wet region (17.3 °C MAT and 603 mm MAP) in the middle of the latitudinal gradient; and a warm-dry region (17.0 °C MAT and 510 mm MAP) in the south. The soils of the farms share a common development from granites, shales and sandy tertiary sediments, are acidic and poor in organic matter, but cover a wide texture gradient (Fig. 1d). In these farms, native pastures are often combined with scattered trees, as is common in Mediterranean and semi-arid rangelands (den Herder et al., 2017; Soliveres et al., 2014). The tree layer is dominated by holm oaks (Quercus ilex L.) with scattered cork oaks (Quercus suber L.) or gall oaks (Quercus faginea Lam.). The herbaceous layer is composed of species typical of Mediterranean pastures and presents a high diversity and proportion of annual C₃ plants (Table S1 in the Supplement). The growing season of the pasture is very limited by the Mediterranean summer drought with the annual species germinating in autumn (early-to mid-October), reaching their peak productivity in mid-spring (late April of the following year) and senescing in June.

Figure 1 Iberian “dehesa” cover according to the CORINE 2018 land cover survey and geographical location of the studied farms (a). Mean annual temperature (b) and precipitation (c) during the period 1980–2018 in the study area. Soil texture in all sampled plots (d).

Dehesa farms are typically managed by extensive continuous grazing, where livestock (mainly cattle) graze freely on large areas following a loosely defined grazing plan. This management has been traditionally practiced in all the farms studied, until, in recent decades, some paddocks were converted to other management practices. As a result, on each farm we selected five paddocks, each with one of the following management regimes:

• Abandoned (Ab): paddocks devoid of grazing for at least the last 10–20 years.

• Continuous grazing (Ct): paddocks where livestock stand most of the year, fed by grazing and supplemented. This is the control treatment in our study as it is the most widespread management in Iberian dehesas.

• Rotational grazing (Ro): paddocks intensively grazed in short periods and with resting periods (i.e. without livestock) lasting for more than 6 months. Rotational grazing has been applied in these paddocks for the last 10–15 years.

• Recent legume sowing (Lr): paddocks where pastures have recently (≤5 years) been sown with legume mixtures. In dehesa farms, legume sowing consists of sowing a mixture of seeds (at a rate of 20 kg ha⁻¹) from various species of annual legumes (pre-inoculated with Rhizobium) such as Trifolium subterraneum L., T. incarnatum L., T. michelianum Savi., T. resupinatum L., T. mutabile Port. and Ornithopus sativus Brot., along with some highly productive annual grass species such as Lolium multiflorum Lam. and Lolium rigidum Gaud., all of which have relatively shallow root systems (Teixeira et al., 2015). These sowings are preceded by surface tillage and phosphorus application to meet the needs of the legumes and stimulate N-fixation (Jongen et al., 2019). Farmers sow legumes mixture only once, because its effect persists over the years thanks to the natural seeds production and the resulting soil seed bank of the sown species. None of the sown plots studied was resown.

• Old legume sowing (Lo): paddocks sown with legume mixtures more than 10 years ago.

The distinction between recent and old legumes sown paddocks allows us to compare short and long-term effects after sowing. In two out of the nine farms, we added an additional control (continuous grazing) paddock close to the legume sowing paddocks because they were at a considerable distance from the other control paddock, limiting their comparability. Additional information on paddocks characteristics is provided in Table S2.

In each paddock, we installed four permanent monitoring plots, in open grassland and away from the direct influence of the trees (at least 10 m away from the trees). This setup allowed us to have a total of 188 sampling plots where we measured soil properties, SOC fractions, and soil microbial communities in 2021 and vegetation traits for 3 years (2021 to 2023).

2.2 Vegetation sampling

We assessed above-ground (ABG) pasture productivity over three growing periods at all plots using $\mathrm{50}\times \mathrm{50}\times \mathrm{50}$ cm grazing exclusion cages. In the 2020–2021 period (hereafter referred to as 2021) we installed one cage on each monitoring plot at the beginning of spring (late March 2021). In the 2021–2022 and 2022–2023 periods (hereafter referred to as 2022 and 2023), we installed the cages at the beginning of the green-up in early November. Standing biomass at the time of cage establishment (t0) was determined by hand clipping within two 25×25 cm squares, randomly placed in the proximity of the cage. In the three years we harvested ABG biomass inside the cages at the beginning of summer in early June (t1), when most of the vegetation was already senescent. Biomass collected both in t0 and t1 was oven-dried at 60 °C for 48 h and weighed. We calculated ABG productivity as the difference between t0 and t1 biomass. Therefore, only spring productivity was measured in 2021, which represents the largest proportion of annual productivity in this system, while a closer estimate of annual productivity was measured in 2022 and 2023.

Plant traits were collected in early May 2021, at the peak of the growing season. We considered a circular area of three meters in diameter around each exclusion cage and within this area, we identified the species present at 25 regularly spaced points using the line intercept method (Godínez-Alvarez et al., 2009), with three concentric circular transects going around the cage. After completing the species inventory for all the sampling plots in a paddock, we collected at least three to ten full individuals, in that same paddock, of each of the identified species, as proposed in the abundance-weighted trait sampling scheme (Carmona et al., 2015). Thus, at least ten individuals of the most abundant species were collected in the same paddock, five individuals for the medium abundant species and three individuals for the rare species. The leaf area (LA), specific leaf area (SLA), specific root length (SLR), leaf dry matter content (LDMC), root dry matter content (RDMC), leaf nitrogen content (LNC) and plant maximum height (see Fig. 2 for more detail) of all the collected individuals was measured following the standard protocols (Pérez-Harguindeguy et al., 2016). In total, more than 10 000 individuals were measured across all the plots. This extensive sampling allowed us to account not only for interspecific differences in trait values, but also for intraspecific variability in trait values across managements, farms and regions. Intraspecific variability has proven to be very important for accurately defining the functional composition of plant communities and their relationship to ecosystem processes (Siefert et al., 2015; Westerband et al., 2021). We also conducted species inventories in the spring of 2022 and 2023, but on these later samplings we only collected and measured the traits of individuals from the species not found in 2021. We used the trait values measured in 2021 for each species in each paddock to impute trait values in 2022 and 2023. The proportion of legumes in each plot was quantified from the species inventories as a measure of the number of plants with N fixation capacity in the communities.

Figure 2 Graphical representation and description of plant functional traits measured. References embedded in the figure are provided below: ¹ Gaudet and Keddy (1988), ² Thomson et al. (2011), ³ Díaz et al. (2016), ⁴ Poorter et al. (2009), ⁵ Kazakou et al. (2009), ⁶ Wilson et al. (1999), ⁷ Freschet et al. (2012), ⁸ Wright et al. (2004), ⁹ Guyonnet et al. (2018), ¹⁰ Roumet et al. (2016), ¹¹ Bergmann et al. (2020), ¹² Kramer-Walter et al. (2016).

Pasture chemical composition was assessed in the three years, after the species inventories and traits samplings. We collected standing biomass by hand-clipping in three 25×25 cm quadrats randomly placed in the monitoring plot, outside the exclusion cages. These samples were dried at 60 °C for 48 h, then grounded and analyzed with Dumas method in DUMATHERM^® N Pro analyzer (C. Gerhardt GmbH & Co. Germany) to obtain the nitrogen content of each sample. Acid detergent lignin (ADL) and fiber (ADF) and neutral detergent fiber (NDF) content were also measured. NDF, ADF and ADL were determined using a fiber analyzer (ANKOM A2000, ANKOM Technology, USA), following the official procedures (Latimer, 2023).

2.3 Soil sampling and soil characteristics measurement

In spring 2021, four soil cores were collected with a push sampler (5 cm diameter) at a depth of 8 cm around each exclusion cage, after removal of above-ground vegetation and litter from the sampled surfaces. The four soil cores were combined to obtain a composite sample from each plot, and an aliquot of 40 g of this sample were sieved (<2 mm), reserved and stored at 4 °C for microbial community analysis in the following days. The remainder of the composite soil samples were air-dried and sieved (<2 mm). Coarse mineral material greater than 2 mm (hereinafter referred to as gravel) was weighed and used for bulk density correction (Eq. 1).

Soil texture (sand, clay and silt content) was determined using a laser diffraction particle size analyzer (Mastersizer 2000, Malvern Instruments Ltd. UK) after dispersion with sodium hexametaphosphate. Soil pH was measured with a glass electrode (soil : water 1 : 2.5) pH meter (CRISON Basic20, Alella, Spain). Soils were extracted with 1M KCL and determined colorimetrically in a Bran + Luebbe Autoanalyzer 3 (Norderstedt, Germany), following the manufacturer's protocol, to measure nitrate (NO${}_{\mathrm{3}}^{-}$) and ammonium (NH${}_{\mathrm{4}}^{+}$) concentrations. Available P (Olsen P) was determined colorimetrically in an Agilent Cary 60 UV-Vis spectrophotometer (Agilent Technologies, USA), after extraction with 0.5 M NaHCO₃ at pH 8.5. Available K, Ca and Mg were determined by ICP-OES (Varian Mod. 720, Palo Alto, California, USA), after extraction with 1 M ammonium acetate (pH 7).

2.4 Soil organic matter fractionation

Representative subsamples of each 2 mm sieved air-dry composite soil sample were used to measure both bulk SOC and total nitrogen (TN) as well as their distribution across distinct physical fractions. First, 10 g aliquots of all samples were ground in a ball mill and the Bernard's calcimeter method (Sherrod et al., 2002) was used to test for the presence of inorganic C. Only four samples (0.02 % of the database) contained traces of carbonates, with CaCO₃ contents between 0.2 % and 0.8 %. These samples were excluded from subsequent analyses. For bulk soil analyses, aliquots of 10g of soil were ground in a ball mill and then analyzed in the DUMATHERM^® N Pro analyzer (C. Gerhardt GmbH & Co. Germany) to determine % C (SOC) and % N (TN). For soil fractionation, we followed a combined size and density procedure as described in Leuthold et al. (2024). Briefly, an aliquot of 6 g of soil was mixed with sodium polytungstate (1.85 g cm⁻³) and shaken reciprocally for 18 h to disperse the aggregates. After dispersion samples were centrifuged for density fractionation and the light particulate organic matter “POM” (<1.85 g cm⁻³) was aspirated from the rest of the soil. We then thoroughly rinsed the residual heavy fraction and separated it by wet sieving into coarse heavy mineral-associated organic matter “chaOM” (>53 µm) and fine mineral-associated organic matter “MAOM” (<53 µm). All fractions were analyzed for % C and % N on an elemental analyzer as described above. Fractionation was accepted when mass recovery was within ±5 %, and C recovery was within ±15 % for samples that did not satisfy one of these conditions, fractionation was repeated.

Figure 3 Violin plot and boxplot (with median and quartiles) and mean values (expressed with labels and red dots) for the ratio of carbon to nitrogen content (C $/$ N ratio) in the soil organic matter fractions (a), and relation between carbon content in the soil organic matter fractions (in g of C per kg of soil) and the total soil organic carbon (SOC) content (b). “chaOM” refers to the coarse heavy mineral-associated organic matter, “lPOM” to the light particulate organic matter and “MAOM” to the fine mineral-associated organic matter. In panel (b), a linear regression line depicts the relationship between IPOM and total SOC, while logistic curves illustrate the relationships of chaOM and MAOM with total SOC.

Since the chaOM and the MAOM shared similar C $/$ N ratios, which were lower than POM C $/$ N ratios (Fig. 3a), we merged these two mineral associated OM fractions and present them together as MAOM (Santos et al., 2024; Zhang et al., 2021). Thus, in this work, POM and MAOM represented the light (<1.85 g cm⁻³) and heavy (>1.85 g cm⁻³) fractions respectively. This contrasts with other studies based on size fractionation, in which light POM and chaOM are pooled together (Cotrufo et al., 2019; Dai et al., 2025; Díaz-Martínez et al., 2024). However, our decision is supported by recent findings of Leuthold et al. (2024) who observed greater chemical similarity between chaOM and MAOM than between chaOM and POM. The pooling of the two mineral associated OM fractions is not expected to modify the results or interpretation of this work, as the C content of chaOM only accounts for an average of 8 % of the total SOC, and its relative importance is similar across all managements (Fig. S1 in the Supplement). Carbon data are presented as SOC, POC and MAOC stocks (Mg ha⁻¹), calculated following (Poeplau et al., 2017):

$$\begin{array}{}\text{(1)}& {\mathrm{OC}}_{\mathrm{stock}}={\mathrm{OC}}_{\mathrm{content}}\times {\displaystyle \frac{{\mathrm{mass}}_{\mathrm{fine}\mathrm{soil}(<\mathrm{2}\mathrm{mm})}}{{\mathrm{Volume}}_{\mathrm{sample}}}}\times {\mathrm{depth}}_{\mathrm{sample}}\end{array}$$

where OC_content is the organic C content (as a proportion) of the soil fraction, mass_{fine soil(< 2 mm)} is the dry weight of the soil excluding gravel and large roots, and considering the dry mass of the aliquot reserved for microbial community analysis. Volume_sample and depth_sample are respectively the volume and depth of the soil sampled with the core.

Additionally, we calculated the proportion of MAOC in total SOC (MAOC $/$ SOC) following Cotrufo et al. (2019)

2.5 Soil microbial communities

Soil microbial communities were characterized by phospholipid fatty acids analysis (PLFA; Willers et al., 2015). A 2 g fine soil aliquot of lyophilized soil was used for lipid extraction with a one-phase chloroform–methanol-phosphate buffer solvent. Phospholipids were separated from non-polar lipids and converted to fatty acid methyl esters (FAMEs), which were then separated by gas chromatography as described in details in Moreno et al. (2021). The sum of all individual PLFAs was used as a proxy for microbial biomass (in nmol PLFAs g⁻¹ of soil). Further, we estimated microbial biomass stocks (in mol ha⁻¹) by substituting the OC_content by the total PLFAs concentration in Eq. (1). Specific PLFAs were used as biomarkers to estimate the relative abundance of broad taxonomic microbial groups, according to their characteristic fatty acids: eukaryote, Gram− and Gram+ bacteria, saprophytic fungi and arbuscular mycorrhizal fungi (Frostegård and Bååth, 1996). The ratios among Gram+ and Gram− bacteria (Gram+ $/$ Gram−) and fungi and bacteria (Fungi $/$ Bacteria) were also calculated to describe the composition of the microbial community.

2.6 Mineral-associated carbon capacity and saturation

To evaluate the degree of MAOC saturation, we used the boundary line approach reported by Georgiou et al. (2022) adjusted with global soil samples as a reference of the maximum observed mineralogical capacity (sensu Georgiou et al., 2025) of our samples (Fig. 4). Therefore, the saturation deficit in each sample was calculated as one minus the ratio of the current C content and the observed maximum C content according to the mineralogical capacity (% clay + silt).

Figure 4 Representation of the 4 main axis of variation in the principal component analysis (PCA) summarizing the vegetation characteristics variables. Panels (a)–(d) show the correlation between the different variables included in the PCA and the new axis. Panels (e)–(h) represent the plant communities characteristics at the extreme of each axis. Representative species of each axis are represented in Fig. S2.

2.7 Data analysis

All analyses were carried out in R 4.3.3 (R Core Team, 2024). Species abundance data and species trait values were used to compute, for each year and plot, the community weighted means (CWM) of all traits as well as the RaoQ multi-trait functional diversity index (Ricotta and Moretti, 2011). In addition, we calculated the species richness as a measure of the taxonomic diversity in each plot. These analyses were performed using the package “FD” (Laliberté and Legendre, 2010). Regarding pasture chemical composition, the hemicellulose content was calculated as NDF minus ADF and the cellulose content as ADF minus ADL, with ADL being considered a good proxy for the lignin content (Van Soest et al., 1991).

For the vegetation analysis, measurements from the three years were combined to obtain mean estimates of each vegetation attribute. For each year, the ABG production and chemical attributes values, as well as the CWM of all the functional traits, were centered and scaled between −1 and 1, so that values closer to −1 represent plots with traits values lower than the mean, while values closer to 1 represent plots above the mean. The three years centered and scaled values were then averaged to obtain estimates of each vegetation characteristic across the study period, as shown in this Eq. (2):

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& \mathrm{Vegetation}{\mathrm{attribute}}_{\mathrm{mean}}=\\ & {\displaystyle \frac{\sum _{\mathrm{year}i=\mathrm{2021}}^{\mathrm{2023}}\mathrm{Vegetation}\mathrm{attribute}(\mathrm{centered}\mathrm{and}\mathrm{scaled}{)}_{\mathrm{year}i}}{\mathrm{3}}}.\end{array}\end{array}$$

All mean values of vegetation characteristics obtained with Eq. (2) were used to build a Principal Component Analysis (PCA). From this PCA, we extracted 4 main axes of variation of the vegetation characteristics with eigenvalues greater than one (Fig. 4).

Nutrients stocks in the soil samples were calculated by substituting the OC_content, in Eq. (1), for the concentration of each nutrient. PCA was used to summarize these nutrient stocks into a single composite index, as all nutrients were positively correlated. The first component of this PCA, which absorbed 46.9 % of the variance, was extracted as a new variable, henceforth called “soil fertility” (Fig. S3). This soil fertility index increases as ammonium, nitrate, P, Ca, Mg, and K stocks in soil increase.

To analyze the direct and indirect effects of management, vegetation characteristics, microbial communities and environmental variables on POC, MAOC and SOC stocks and MAOC $/$ SOC ratio we built a structural equation model (SEM) using the “piecewiseSEM” package (Lefcheck, 2016). These models assess the extent to which a defined structure of causal relationships between variables fits the actual correlations in the data. In this causal structure, direct effects describe the influence of one variable on another without intermediaries, whereas indirect effects occur when a variable affects another through its impact on a third variable. We started from the assumption that management, climate and soil properties could affect SOC, POC and MAOC stocks and the MAOC $/$ SOC directly or indirectly, by modifying microbial composition and vegetation characteristics. This provided us with a theoretical basis for the construction of the SEMs. The linear sub-models within these SEMs were fitted with the “lme4” package (Bates et al., 2015) using the farms as random factors to avoid spatial pseudo-replication. Model fitting was performed according to Zuur et al. (2009) selecting the model with the lowest corrected Akaike Information Criterion (AICc) value after adjusting for random and fixed factors. All predictor variables used in the model selection are summarised in Table 1. While gravel content is correlated to BD, its inclusion as predictor in the models helped to control for inter-site variability in gravel content. Interactions between management and climate or soil texture were also included during model selection. All model assumptions were checked and satisfactorily met. Predictor variables were scaled and centered prior to inclusion in the models.

Table 1 List of the variables and their units, which were included as fixed factors during the structural equation model (SEM) sub-models selection.

To calculate the total effect of each explanatory factor (i.e. the sum of direct and indirect effects) on SOC, POC and MAOC stocks, and MAOC $/$ SOC ratio in the SEMs, we used the “semEff” package (Murphy, 2020). This package generates estimates and confidence intervals of the total effects of the explanatory factors across multiple permutations of the data (9999 bootstraps in our case).

Given that average bulk density values partly differed among managements (Fig. S4), and that this variation may affect carbon stocks, the same analysis procedure used for POC, MAOC and SOC stocks was used to analyze the change in POC, MAOC and SOC contents.

3 Results

The mean values of POC, MAOC and SOC stocks were 9.8±4.3, 7.1±3.1 and 16.5±6.3 Mg C ha⁻¹ respectively. On average, MAOC represented 44 % of the total SOC. The MAOC content of all the soils analyzed was well below the theoretical saturation level established by Georgiou et al. (2022) in relation to the clay + silt content (Fig. 5). Mean MAOC saturation deficit was 80 %. Soil organic carbon fractions and stocks were affected directly and indirectly by mean plant traits, microbial communities, soil properties and management, and indirectly by climatic conditions (Fig. 6). Plant diversity indices and interactions between management and climate or soil texture were not retained in the model selection due to lack of significance.

Figure 5 Relation between the mineral-associated organic carbon (MAOC) and the percent Clay + Silt (CS) in the soil samples and the boundary line (in blue) with confidence intervals adjusted by Georgiou et al. (2022) for high-activity mineral soils, indicating the maximum observed mineralogical capacity for each CS content. The equation for the boundary line is provided in the panel. Points are colored based on their MAOC deficit (1 − MAOC content/MAOC maximum observed capacity).

Figure 6 Structural equation model representation. Factors included in the model are grouped by factor type (management, climate, vegetation traits, soil properties, microbial communities). Arrows between groups of factors indicate significant relationships between any of the factors included in both groups. The width of these arrows is proportional to the mean absolute size of the estimates between the factors in the groups. The plots embedded into these arrows show the standardized estimates of the significant relationships between the factors connected by the arrow. Negative standardized estimates are represented in red, and positive ones in blue. The size of the estimate circles represents the absolute value of the standardized estimate. Causal relationships between factors in the same group are represented by individual arrows.

Total microbial biomass, estimated as total PLFA stocks, averaged 66.3±25.9 mol ha⁻¹, with mean values for Fungi $/$ Bacteria and Gram+ $/$ Gram− ratios of 0.08±0.03 and 1.1±0.2, respectively, although there were some differences between treatments (Fig. S5). Microbial biomass was positively correlated with POC, MAOC and SOC stocks (Figs. 6f and 7a–c). The Fungi $/$ Bacteria ratio was negatively correlated with MAOC stocks (Figs. 6f and 7b) meaning that these stocks were higher in soils with a greater bacterial predominance. The Gram+ $/$ Gram− ratio had a significant direct negative effect on POC, MAOC and SOC stocks (Figs. 6f and 7a–c).

Figure 7 Direct, indirect and total standardized effects of all studied variables included in the structural equation model (Fig. 6) over the (a) particulate organic carbon (POC), (b) mineral-associated organic carbon (MAOC), and (c) soil organic carbon (SOC) stocks and (d) the relative MAOC abundance (MAOC $/$ SOC). Bars indicate direct (dark colors) and indirect (light colors) effects, and the black points-ranges indicate the total (i.e. direct + indirect) effect (with its 95 % confidence interval). Stars over the total effect values indicate significant effects at a level of 0.05.

The first axis of the PCA of the vegetation characteristics, which reflects the N-C stoichiometry of the vegetation, had a negative direct effect on POC, MAOC and SOC stocks (Figs. 6h and 7a–c). The second axis, related to the plant economic spectrum (PES), had a positive direct effect on SOC stocks (Figs. 6h and 7c). In addition, the PES axis was positively correlated with the microbial biomass, thus having a positive indirect effect, and a significant total effect on POC, MAOC and SOC stocks (Figs. 6e and 7a–c). The third axis of the PCA (P-L axis), positively correlated with the lignin content and negatively related to vegetation productivity, was negatively correlated with the POC, MAOC and SOC stocks and the microbial biomass (Figs. 6e, h; and 7a–c). The fourth axis, reporting the plant size, was positively correlated with the MAOC stocks, the MAOC $/$ SOC ratio and the Fungi $/$ Bacteria ratio (Figs. 6e, h; and 7a, d). No index of plant functional and taxonomic diversity was retained during model selection. On the other hand, soil fertility had a direct and indirect positive effect on POC, MAOC and SOC stocks (Fig. 6g). Soil fertility was negatively correlated with the Gram+ $/$ Gram− ratio, the P-L axis and the N-C axis, and positively correlated with the PES axis of the vegetation PCA (Fig. 6i and j).

Plot under continuous grazing had a mean POC, MAOC and SOC stocks of 9.1, 6.8 and 16.2 Mg C ha⁻¹ respectively, and a MAOC $/$ SOC ratio of 0.44 (Fig. 8). Rotational grazing significantly increased soil fertility compared to continuous grazing (Fig. 6k). This led to a significant indirect effect of rotational grazing on POC, MAOC and SOC stocks (Fig. 7b and c). POC, MAOC and SOC stocks under rotational grazing had a mean value of 10.5, 7.8 and 18.5 Mg C ha⁻¹, respectively, which were 15 % (in the case of POC and MAOC) and 14 % (for SOC stocks) higher than mean values in continuous grazing (Fig. 8b and c). Recent legume sowing had a negative direct effect on POC, MAOC and SOC (Fig. 6b). However, Lr also increased significantly the plant size axis, the microbial biomass and the soil fertility and decreased the plant N-C axis compared to Ct (Fig. 6a, c and k). These changes resulted in a positive indirect effect of Lr over POC, MAOC and SOC stocks and a null total effect (Fig. 7a–d). Grazing abandonment increased the vegetation N-C axis and reduced the microbial biomass compared to Ct (Fig. 6a and c), resulting in a negative indirect effect on POC, MAOC and SOC stocks and a significant negative total effect on POC stocks (Fig. 7a–c). Thus, the mean POC stock on abandoned plots was 8.1 Mg C ha⁻¹, 11 % lower than in Ct plots (Fig. 8a).

Figure 8 Violin plot and boxplot (with median and quartiles) for the (a) particulate organic carbon (POC), (b) mineral-associated organic carbon (MAOC), and (c) soil organic carbon (SOC) stocks and (d) the relative MAOC abundance (MAOC $/$ SOC) in each management (Ct = continuous grazing; Ab = grazing abandonment; Lr = recent legume sowing; Lo = old legume sowing; Ro = rotational grazing). The red point and the label indicate mean values predicted by the structure equation model (Fig. 6). Lower case letters indicate significant differences (p<0.05) between managements.

Sand content in soil was negatively correlated with MAOC stocks and MAOC $/$ SOC ratio (Fig. 6j). Sand content was also negatively correlated with soil fertility and the PES axis and positively correlated with the P-L and N-C axis (Fig. 6i), thus having a negative indirect effect on POC, MAOC and SOC stocks (Fig. 7a–c). As expected, due to its negative correlation with BD, gravel content was also negatively correlated with microbial biomass (Fig. 6j), leading to a significant negative indirect effect on POC, MAOC and SOC stocks (Fig. 7a–c). MAT was positively correlated with the N-C axis and negatively correlated with the soil fertility index (Fig. 6d and l). This implied a negative indirect effect of MAT over POC, MAOC and SOC stocks (Fig. 7a–c). On the other hand, MAP was negatively correlated with the N-C and the size axis (Fig. 6d). Thus, MAP had a positive indirect effect on POC, MAOC and SOC stocks (Fig. 7a–c and e).

4 Discussion

Our results highlight the potential of management to control carbon storage in semi-arid grasslands. In this sense, rotational grazing arises as a promising tool for enhancing long-term carbon storage in the topsoil. We also identified several pathways by which management influences SOC, POC and MAOC stocks, showing the importance of changes in vegetation and microbial composition. Furthermore, management effects were consistent across the wide soil-climate gradient examined, enabling generalization of our results to a broad set of semi-arid grasslands in various environmental contexts. It should be noted that the predictive capacity of POC, MAOC and SOC stocks models was high (62 %, 69 % 73 % of the variance explained by the fixed factors, respectively), indicating the robustness of the results presented and their importance for expanding the understanding of soil carbon dynamics in grasslands. Our results are limited to the upper topsoil (the first 8 cm), and it would be necessary to analyze deeper layers to fully understand the processes of SOC formation and stabilization and the effects of management changes. However this top layer is the most important for mediterranean grassland functioning as it contains the majority of roots and therefore most of the microbial, nutrients and water dynamics (Acosta-Gallo et al., 2011; Moreno et al., 2005). Furthermore, changes in management primarily affect the SOC of the topsoil layer, especially in the short term (Ward et al., 2016).

4.1 Soil organic carbon fractions, stocks, and saturation levels

SOC stocks and contents (Figs. 8 and S4) in the paddocks studied (16.5 Mg C ha⁻¹ and 19.5 g C kg⁻¹ on average, respectively) were consistent with values found in the first 7–10 cm of soil in other grasslands under similar climatic conditions (Díaz-Martínez et al., 2024; Oggioni et al., 2020; Parras-Alcántara et al., 2014). On the other hand, the MAOC $/$ SOC ratio of our soils (0.44 on average) was lower than the values generally reported in both global and semi-arid grasslands (Hansen et al., 2024; Rocci et al., 2022; Shi et al., 2024), which are around 0.6–0.7. The shallow sampling depth in this work (8 cm) may influence the observed MAOC $/$ SOC ratios, as POM decrease more sharply than MAOM with soil depth (Galluzzi et al., 2025; Sanderman et al., 2021; Tang et al., 2024). For instance, Plaza et al. (2022) reported MAOC $/$ SOC ratios between 0.3 and 0.4 in the first 10 cm of soil, while Cappai et al. (2017) found MAOC $/$ SOC ratios around 0.7 in the first 20 cm, both in Mediterranean grasslands.

Carbon concentrations in the fine soil fraction (clay + silt) of the studied grasslands were well below the maximum C capacity observed in previous studies (Cotrufo et al., 2019; Georgiou et al., 2022). However, we detected a certain limit to MAOC accumulation, as its content remained stable above SOC contents of 30 g kg⁻¹ and stayed below 20 g kg⁻¹ even when SOC content exceeded 60 g kg⁻¹, following a saturation curve (Fig. 3b). This result suggests that MAOC accrual in this system is mainly limited by environmental conditions rather than by the mineralogical capacity of the soil or the amount of inputs. In this sense, the maximum C content in the clay + silt fraction observed in this work (around 17 g C kg⁻¹ clay + silt) likely represent the effective capacity of the ecosystem (sensu Georgiou et al., 2025) to stabilize C under current environmental conditions. Although this constraint may limit management opportunities to increase MAOC stocks through management, most plots had MAOC contents well below this effective capacity, indicating substantial room for improvement. Indeed, we observed significant management effects on the MAOC fraction. Importantly, management decisions should not rely solely on MAOC accumulation, as POC can also persist for long periods of time, acts as a precursor of MAOC, and has no known upper limit for accumulation (Angst et al., 2023).

4.2 Soil microbial communities regulate SOC storage

Total microbial biomass had a similar and substantial positive effect on POC, MAOC and total SOC stocks, such that the MAOC $/$ SOC ratio was unaffected. Typically, microbial biomass C only represents about 2 % of SOC (Xu et al., 2013; Yao et al., 2000), but it is closely correlated with the accumulation of microbial necromass C (Hou et al., 2024). The latter can account for more than 50 % of SOC (Liang et al., 2019), with a similar contribution to the POC and MAOC fractions (Zhang et al., 2024). The positive effect of microbial biomass on both POC and MAOC, as well as the relatively low C : N ratio of these fractions found in this study (Fig. 3a) compared to other works (Chang et al., 2024; Yu et al., 2022), indicate a prominent role of microbial transformation of plant inputs in SOC formation in semi-arid grasslands, which would also explain the relatively low SOC stocks of these ecosystems.

Fungi $/$ bacteria ratio was negatively correlated with MAOC fraction, as fungi residues tend to contribute more to POC (Griepentrog et al., 2014; Lavallee et al., 2020; Tang et al., 2023), but no significant effect of fungi $/$ bacteria ratio over POC stocks or MAOC $/$ SOC ratio was observed. Gram+ $/$ Gram− ratio was negatively correlated with POC, MAOC and SOC stock, as observed in previous studies (Khatri-Chhetri et al., 2024). Gram− bacteria are more dependent on plant C inputs, whereas Gram+ bacteria tend to use more of the organic C already present in the soil (Fanin et al., 2019; Kramer and Gleixner, 2008; Waldrop and Firestone, 2004). Thus, the proliferation of Gram+ bacteria may promote the decomposition of pre-existing SOM (Klumpp et al., 2009).

4.3 Plant–soil interactions and their effect on soil carbon stocks

The analysis of vegetation functional traits and chemical composition revealed four main axes of variation. Two of these axes (axes 1 and 4) correspond to the spectrum of plant form and function, found in several global vegetation analyses (Díaz et al., 2016; Weigelt et al., 2021). This two-dimensional spectrum is defined by the leaf and root economic gradient (Kramer-Walter et al., 2016; Wright et al., 2004), that moves from slow-growing and resource-conserving to fast-growing and resource-acquisitive species (Wright et al., 2004), and the size gradient, which increases with increasing plant height and leaf area and reflects the competitive ability of plants for the light (Díaz et al., 2016).The other two axes (axes 2 and 3) are more related to the aboveground chemical composition of the vegetation. Axis 1 reflects the N-C stoichiometry of plant tissues, driven mainly by the ABG C $/$ N ratio, both of which are obviously negatively correlated with certain traits such as LNC and the proportion of legumes. The N-fixing capacity of legumes has a fertilising effect on companion plants, increasing its tissues nitrogen content (Pino et al., 2016). In addition, legumes themselves tend to have a high concentration of nitrogen in leaves and stems (Carranca et al., 2015). Axis 3 is positively related to the lignin content of the vegetation. This independence between lignin and nitrogen content in plant litter has already been observed (Cornwell et al., 2008). ABG productivity of vegetation was correlated with most of these axes, especially with axis 3, being higher in more acquisitive, less lignified and bigger vegetation, as found in other works (Laliberté and Tylianakis, 2012; Zhang et al., 2017).

POC, MAOC and SOC stocks were higher in communities dominated by resource-acquisitive, highly productive plants with low lignin content and low C $/$ N ratio. These communities are expected to provide a higher level of plant inputs, thereby increasing the incorporation of organic matter into the soil (King et al., 2023; Mortensen et al., 2025; Zhou et al., 2024). Increased inputs may come from higher plant litter production but also from increased root exudates, which tend to be higher in acquisitive plants (Guyonnet et al., 2018). Part of the effects of plant traits on SOC was mediated by increases in microbial biomass, which would also benefit from higher levels of plant inputs and increased root exudation (Eisenhauer et al., 2017). In addition, inputs with low lignin and C $/$ N ratios are degraded more efficiently, reducing C losses and promoting long-term SOC storage (Cotrufo and Lavallee, 2022; Ridgeway et al., 2022). However, based on previous work, we would have expected a direct contribution of plant structural input to POC (Cotrufo et al., 2022), and thus a higher proportion of POC in communities that produce more recalcitrant litter (Cheng et al., 2023; Haddix et al., 2016; Mortensen et al., 2025). This was not the case in this work, where the MAOC $/$ SOC ratio remained unaffected by the chemical composition of the ABG vegetation tissues. This unexpected result could be explained by the fact that in these systems, POC also appears to be the product of microbial transformation of plant inputs, as suggested by its positive relationship with microbial biomass and its relatively low C $/$ N ratio. In addition, photodegradation, an important driver of litter degradation in semi-arid ecosystems (Austin and Vivanco, 2006), can promote litter lignin biodegradability and the production of litter soluble compounds that are readily accessible to soil microbes (Wang et al., 2015), thus reducing the influence of vegetation chemical properties on POC and MAOC formation. Grazing also decouples the quality of plant tissues and the final quality of inputs to the soil, as livestock excreta have a lower C $/$ N ratio than the plant material consumed (Soussana and Lemaire, 2014). We also observed that taller and larger plants (high values on the size axis) were associated with higher MAOC stocks and MAOC $/$ SOC ratios, although the mechanism driving this correlation was unclear. Generally, plant height is positively correlated with shoot : root ratio (Li et al., 2008), and several studies have found a higher contribution of shoots, rather than roots, to the MAOC fraction due to the higher recalcitrance of root tissues (Huang et al., 2021; Lavallee et al., 2018; Ridgeway et al., 2022). However, rhizodeposition, which is closely linked to root growth and turnover, has been shown to promote MAOC formation over POC (Berenstecher et al., 2023; Engedal et al., 2023; Villarino et al., 2021). On the other hand, a greater accumulation of standing litter, rather than surface litter, might be expected in communities with bigger plants, and some studies in semi-arid grasslands have observed higher rates of microbial degradation and release of soluble compounds (thus contributing more to MAOC) in standing litter, compared to surface litter, due to greater retention of night-time moisture in the former (Gliksman et al., 2018; Wang et al., 2017). In any case, our results and their interpretation are limited by the correlation in our vegetation data between input quantity (ABG productivity) and input quality, which could mask the effect of the latter on SOC formation.

No significant effects of taxonomic or functional diversity on POC, MAOC or SOC stocks were found, which could be explained by the negative correlation of these biodiversity indicators with the ABG productivity and soil fertility in our plots (Fig. S7). This negative correlation has been observed in other natural grassland research and may be due to the dominance of some highly productive species in nutrient-rich environments (Feßel et al., 2016; Helsen et al., 2014; Luo et al., 2019; Rolo et al., 2016). In this sense, our results would be consistent with previous research that has found a greater importance of functional identity (i.e. mean traits values) relative to functional diversity in predicting ecosystem processes (Mokany et al., 2008; Zhang et al., 2017).

Soil fertility plays a central role in our models, increasing POC, MAOC and SOC stocks. Previous studies in grasslands have already observed a positive correlation between POC and MAOC and soil nutrients content (Wang et al., 2025). Global N addition experiments have also showed to increase topsoil C storage in grasslands (Hu et al., 2024; Liu et al., 2023). In this work, we observe a direct effect of soil fertility on C storage, but we also identify some indirect effects mediated by changes in vegetation and microbial communities. More resource-acquisitive and productive plant communities, with lignin content and lower C $/$ N ratio, were found at higher soil fertility. Resource-acquisitive plants are expected to be more competitive in nutrient-rich environments (Daou et al., 2021; Ordoñez et al., 2009). The absence of nutrient limitation also promotes plant productivity (Fay et al., 2015). In addition, soil fertility was negatively correlated with the Gram+ $/$ Gram− ratio. Gram- bacteria have a more copiotroph strategy than Gram+ bacteria and therefore benefit more from nutrient-rich environments (De Deyn et al., 2010; Luo et al., 2020; Zhong et al., 2010). The direct positive effect of soil fertility on C storage in soils may be explained by changes in other factors that we did not measure. For instance, nutrient addition have been shown to promote below-ground productivity in Mediterranean woody pastures (Nair et al., 2019) and in global grasslands (Liu et al., 2023; Lu et al., 2011). Moreover, nutrient availability in soils tends to increase CUE (Poeplau et al., 2019; Spohn et al., 2016) and to alleviate the need for microbial N-mining, reducing C losses and old SOM decomposition (Blagodatskaya and Kuzyakov, 2008).

4.4 Management effects

Rotational grazing increased both POC, MAOC and SOC stocks by 15 %, 15 % and 14 %, respectively, compared to conventional continuous grazing. Previous studies found the same positive effect of rotational grazing on MAOC (Khatri-Chhetri et al., 2024; Mosier et al., 2021; Stanley et al., 2025) and total SOC content (Byrnes et al., 2018; Conant et al., 2017; Liu et al., 2024; Phukubye et al., 2022; Teague et al., 2011). According to our results, this effect is mainly explained by the higher soil fertility observed in rotational grazing plots, cascading in the above explained effects on plant traits and microbial transformation leading to higher POC, MAOC and SOC stocks. Soil fertility may be higher under rotational grazing than under continuous grazing because fecal and urine excretions tend to be less spatially clustered in rotationally grazed paddocks, resulting in more homogeneous fertilization of the entire paddock (Augustine et al., 2023; Dubeux et al., 2014). In contrast, in continuous grazing, excreta tend to accumulate in areas of highest animal use, such as near feeders or water points (Tate et al., 2003), which are areas that have been avoided in our sampling. In turn, the homogeneous grazing and fertilization maintains a greater amount of ground covered by vegetation (Stanley et al., 2025; Teague et al., 2004), which limits topsoil losses through erosion (Sanjari et al., 2009).

Surprisingly legume sowings had a negative direct effect on soil C stocks, but a positive indirect effect, resulting in a non-significant total effect. The negative impact was more evident in recent sowings, especially on the MAOC fraction, and could be attributed to the impact of the pre-sowing tillage. Previous works have observed a reduction in SOC stocks after tilling in mediterranean grasslands (Parras-Alcántara et al., 2015; Uribe et al., 2015). This would explain why the direct effect of legume sowing in C storage is more negative in recent compared to old-sowed paddocks. The effect of tillage on SOC stocks may arise from changes in soil bulk density rather than C content, resulting in misleading conclusions (Du et al., 2017; Rovira et al., 2022). However, in this work the analysis of POC, MAOC and SOC contents (in g kg⁻¹ bulk soil; Figs. S6 and S8) shows similar effects of legume sowing as those observed for C stocks (Fig. 7a–c). In the recently sown plots, the negative impact in C stocks was countered by an increase in soil fertility, plant productivity and plant nitrogen content, possibly due to increased legume cover (Gómez-Rey et al., 2012; Gou et al., 2023; Hernández-Esteban et al., 2019). In addition, recent legume sowing also implied an increase in microbial biomass and a decrease in the Gram+ $/$ Gram− ratio, as observed in previous studies (Moreno et al., 2021).

Grazing abandonment led to a 12 % reduction in the POC stocks, compared to continuous grazing. This result contradicts the observations of previous studies in grasslands, which reported positive effects of grazing exclusion on SOC storage, both globally (Eze et al., 2018) and particularly in semi-arid climates (Cheng et al., 2016; Gao et al., 2025; Yu et al., 2021). However McSherry and Ritchie (2013) meta-analysis found that light grazing could promote SOC stocks in C₃ dominated grasslands (such as our study grassland). Other global and regional studies also observed null or even positive effects of light grazing, compared to grazing exclusion, in SOC stocks in grasslands (Liu et al., 2024; Zhou et al., 2017). Our results are in line with previous studies in Mediterranean woody pastures that found a decrease in SOC storage with grazing abandonment (Oggioni et al., 2020; Peco et al., 2017). The fact that grazing abandonment particularly affected the POC fraction is surprising, since in abandoned paddocks there were a greater accumulation of plant structural biomass, that should promote POC accrual (Cotrufo et al., 2022). According to our results, the negative effect of grazing exclusion on POC stocks was mainly mediated by a reduction in microbial biomass, and a proliferation of plants with higher C $/$ N ratios. Once again pointing to the microbially transformed nature of POC in these soils. Grazing has been proven to increase microbial activity and biomass (Bardgett et al., 2001; Zhou et al., 2017), in part by increasing root exudation (Hamilton et al., 2008; Wilson et al., 2018). Previous studies have also observed an increase in vegetation C $/$ N ratio with cessation of grazing (He et al., 2020; Wang et al., 2016), which could be explained by the accumulation of older plant tissues and senescent standing biomass, with higher C $/$ N than green biomass and new leaves (Sanaullah et al., 2010).

Importantly, we observed no interactions between management effects and climate or soil texture variables, showing that the observed management effects were consistent across the entire environmental gradient. Although such interactions have been observed for wider climatic gradients (Byrnes et al., 2018; Phukubye et al., 2022), our results show that management practices such as rotational grazing can have net positive impacts on SOC storage in grasslands with varying conditions within semi-arid regions.

4.5 Climate effects and implications for global warming

Even within the small climatic gradient covered in this study, we found lower stocks of POC, MAOC and SOC in the warmer areas, as observed in previous research (Díaz-Martínez et al., 2024; Georgiou et al., 2024; Hansen et al., 2024; Shi et al., 2024). According to our results, these negative effects were entirely mediated by a decrease in soil fertility and vegetation input quality (increased N-C axis) with increasing temperature. Increased aridity has been shown to be associated with the alteration of several ecosystem processes that regulate nutrient cycling and availability (Berdugo et al., 2020; Moreno-Jiménez et al., 2019). Increased fiber for tissue protection, typical from plants growing at higher temperatures, may explain the increase in vegetation C $/$ N ratios (Arroyo et al., 2024). On the other hand, annual precipitation exerted a small positive effect on POC and SOC stocks, mediated by the reduction of C $/$ N ratios in vegetation in wetter farms. Based on future projections of increasing temperatures and decreasing annual precipitation in our study area (IPCC, 2023b) we could expect significant reductions in SOC stocks in Mediterranean grasslands in the next decades, affecting both the POC and MAOC fractions.

5 Conclusions

Implementing rotational grazing can increase both topsoil POC and MAOC stocks, whereas the abandonment of grazing may reduce the carbon storage capacity of Mediterranean grasslands. We found that most of these management effects were mediated by changes in soil nutrients stocks, microbial communities and vegetation morpho-chemical traits, which were, in turn, the main drivers of changes in SOC stocks and fractions. Specifically, more fertile soils with higher microbial biomass, lower Gram+ $/$ Gram− ratios, and more productive plant communities, with acquisitive traits and nitrogen-rich tissues, promoted the accrual of both POC and MAOC. In this regard, and as a novelty, this work has also proven the usefulness of plant functional traits as tools for the study of plant-soil interactions and SOC formation dynamics. Considering global change trends, our findings suggest a potential loss of soil carbon in the studied grasslands over the coming decades due to grazing abandonment and increasing temperatures. Selecting management approaches that mitigate or counteract these losses is vital to maintaining the fertility, productivity and functioning of semi-arid grasslands.