Soil organic carbon (SOC) represents one of the largest carbon reservoirs in the Earth system and plays a critical role in regulating the global climate, supporting nutrient cycling, and sustaining the productivity of agroecosystems (Lal et al., 2021; Kopittke et al., 2022; Da Gama, 2023). The stability and dynamics of SOC are modulated by biogeochemical processes, which are, in turn, influenced by land use and management practices as well as local environmental conditions (Zhang et al., 2024; Maqbool et al., 2025). Understanding how SOC stocks respond to different management practices is essential for evaluating ecosystem carbon-sink potential and their contribution to mitigating global climate change (Ramesh et al., 2019; Mäkipää et al., 2023; Schreiner-McGraw et al., 2024).

Studies indicate that most of the organic carbon in the soil is associated with macroaggregates, this includes carbon that is trapped within these aggregates, this carbon represents a substantial proportion of the total stock and responds more quickly to changes in management (Howlett et al., 2011; Gelaw et al., 2015). By contrast, carbon associated with finer fractions, microaggregates and silt-clay complex, tends to have lower absolute concentrations, but is more stable due to organomineral protection and a longer residence time in the soil (Mustafa et al., 2021; Delahaie et al., 2024). This structural understanding of organic matter supports the hypothesis that changes in disturbance intensity, such as grazing, primarily affect the most labile and superficially distributed compartments, while the more stable pools remain relatively resilient (Dong et al., 2020; Encarnation et al., 2025). In this context, silvopastoral systems have been described as being structurally resilient to moderate disturbances, as they maintain the integrity of soil aggregates and the biogeochemical functionality of the soil (Mosquera-Losada et al., 2018).

When managed sustainably, grazing can contribute to the maintenance or even enhancement of soil organic carbon (SOC) through several mechanisms, including the incorporation of plant residues into the soil, stimulation of root turnover, and the addition of organic inputs via animal excreta (McDonald et al., 2023; Ritchie, 2020; Gross and Glaser, 2021). However, excessive grazing pressure may lead to soil compaction, reduced water infiltration, and loss of vegetation cover, which can ultimately reduce SOC stocks and alter soil functioning (Abdalla et al., 2018; Donovan and Monaghan, 2021; Centeri, 2022; Fonseca et al., 2023). Consequently, the overall impact of grazing on soil carbon dynamics is highly context-dependent, being influenced by factors such as stocking rate, grazing season, management practices, and local environmental conditions (Eze et al., 2018; Santos et al., 2022; Stanley et al., 2024).

However, several studies have also reported limited or non-significant effects of grazing on soil carbon and nitrogen stocks, particularly in extensive or well-managed silvopastoral systems (Piñeiro et al., 2010; Ritchie, 2020; Abagandura et al., 2024). Under moderate grazing pressure and with the maintenance of vegetation cover, soil disturbance tends to be reduced, allowing soil organic matter stocks to remain relatively stable over time (Mosier et al., 2021; Zhou et al., 2022). In such situations, the balance between organic matter inputs from plant residues and animal excreta and potential losses caused by trampling or biomass removal may result in neutral net effects on soil carbon and nitrogen dynamics (Coonan et al., 2020; Zhu et al., 2021; Li et al., 2025). These findings highlight the importance of evaluating grazing impacts within specific ecological and management contexts, as soil responses may vary considerably depending on site conditions, grazing intensity, and soil characteristics (Eze et al., 2018; Liu et al., 2021).

In the face of climate change and growing pressure on agri-food systems, the integration of grazing practices that align with carbon sequestration is emerging as a promising strategy for reconciling sustainable livestock production and climate regulation (Godde et al., 2020; Rui et al., 2022; Stanley et al., 2024). However, grazing intensity plays an important role. Low stocking rates in silvopastoral systems tend to promote greater carbon sequestration in the upper soil layers compared with intensive grazing or grazing exclusion (Ferreiro-Domínguez et al., 2016). In addition, the integration of trees into pasturelands through silvopastoral management aligns livestock productivity with long-term strategies for soil carbon conservation and climate change mitigation (Torres et al., 2022; Aryal et al., 2022).

In northeastern Portugal, native oak forests such as Quercus rotundifolia and Quercus pyrenaica are renowned for their high ecological and cultural relevance (Carvalho, 2011; Pérez-Luque et al., 2020; Pousa et al., 2025). Hese woodlands support a variety of traditional practices, with extensive grazing playing a significant part in shaping the region’s rural landscape (Guadilla-Sáez et al., 2019; Plieninger et al., 2021; Castillo-Garcia et al., 2022). Maintaining an appropriate balance between conservation goals and pastoral use is essential to preserving ecosystem services, including soil protection, nutrient cycling, biodiversity conservation, and the provision of silvopastoral products (Díaz-Pereira et al., 2020; Lecegui et al., 2022; Amorim et al., 2023). In addition, grazing plays an important socio-economic role in north-eastern Portugal, particularly for small and medium-sized family farms, contributing to income generation, rural population retention and the preservation of traditional knowledge (Torres-Manso et al., 2017; Castro et al., 2021; Rego et al., 2022; Sharifian et al., 2022).

Despite the ecological and climatic importance of the subject, significant gaps remain in the literature concerning the impact of grazing on soil carbon stocks in native forests, particularly in oak forests in Portugal. Most available studies focus on permanent pastures (Deosaran et al., 2024), intensive agricultural systems (Rodrigues et al., 2019), or forest ecosystems evaluated in isolation (Patrício et al., 2023). Integrated systems, such as silvopastoral systems, remain relatively understudied. This hinders our understanding of the mechanisms that regulate soil carbon dynamics and stability in environments where trees, shrubs, herbaceous plants, and pastoral management interact simultaneously.

Therefore, this study aimed to evaluate the effects of grazing on soil carbon and nitrogen dynamics and their distribution among particle-size fractions, in two Mediterranean oak silvopastoral systems in northeastern Portugal, namely Quercus pyrenaica and Quercus rotundifolia.

Principal component analysis (PCA) of soil properties in ungrazed plots showed that the first two components explained 86.8% of the total variance (Dim1: 59.7%; Dim2: 27.1%) (Figure 2). Dim1 represented a clear gradient of soil organic matter accumulation, strongly driven by carbon stocks (C Mg ha⁻¹), carbon content (C g kg⁻¹) and nitrogen content (N g kg⁻¹), which exhibited the highest contributions (indicated by longer vectors and warmer colours). In contrast, bulk density was negatively associated with this axis, indicating that soils with higher compaction were characterised by lower C and N levels. Dim2 was primarily associated with soil moisture, reflecting a secondary environmental gradient independent from carbon and nitrogen dynamics. The orthogonal orientation of moisture relative to C and N variables suggests that water availability was not the main driver of nutrient distribution in these systems. The ordination revealed a clear separation between stand types. Samples under Q. pyrenaica were predominantly associated with the positive side of Dim1, indicating higher soil carbon and nitrogen contents and stocks. In contrast, samples under Q. rotundifolia were grouped towards the negative or central region of the axis, reflecting soils with lower nutrient levels and higher bulk density. This pattern indicates that vegetation type is the main factor structuring soil C and N distribution, overriding potential variability associated with soil moisture.

Following the grazing period, the PCA revealed a consistent ordination structure, with the first two principal components explaining 62.1% and 24.8% of the total variance, respectively (Figure 3). As observed prior to grazing, Dim1 represented a strong gradient of soil carbon and nitrogen status, driven primarily by carbon stock (C Mg ha⁻¹), carbon content (C g kg⁻¹), and nitrogen content (N g kg⁻¹), which exhibited the highest contributions (longer vectors with warmer colours). Bulk density remained negatively associated with this axis, reinforcing the inverse relationship between soil compaction and organic matter accumulation. Dim2 was mainly associated with soil moisture, indicating a secondary environmental gradient largely independent of the distribution of carbon and nitrogen. The near-orthogonal orientation between moisture and the C–N variables suggests that water availability did not directly control the spatial variability of soil organic matter in the studied systems. The distribution of samples after grazing closely mirrored the pattern observed in ungrazed conditions, with a clear separation between stand types along Dim1. Samples under Q. pyrenaica remained associated with higher carbon and nitrogen contents and stocks, whereas those under Q. rotundifolia were predominantly located in regions characterised by lower nutrient levels and higher bulk density. No substantial temporal shifts were detected in the ungrazed sites during the study period, with control samples maintaining a relatively stable position and dispersion in the ordination space compared with the grazed treatments.

Overall, the stability of the ordination pattern between sampling periods indicates that grazing, under the conditions evaluated, did not substantially alter the multivariate structure of soil properties. Instead, vegetation type and inherent soil characteristics remained the dominant factors controlling the distribution of soil carbon and nitrogen.

The sandy fraction showed the highest levels of C and N in the soil. It was found that most of the variables evaluated were significantly influenced by habitat type and soil depth, while grazing had no statistically significant effect (Table 1; Tables 2 and 3).

Soil carbon (g kg⁻¹) was significantly influenced by habitat type (Q. rotundifolia vs. Q. pyrenaica) (p < 0.001) and soil depth (p < 0.001) in all fractions analysed (Table 1). Grazing had no significant effect on soil carbon fractions (p > 0.05). None of the interactions between the factors evaluated were significant for this variable (p > 0.05).

The Tukey test indicated significant differences mainly between habitats, with higher soil carbon contents (C g kg⁻¹) in areas dominated by Q. pyrenaica compared to Q. rotundifolia (Figure 4). In the surface layer (10 cm), C values ranged from approximately 35.00 g kg⁻¹ (±2.07 SE) under Q. pyrenaica to 24.60 g kg⁻¹ (±1.45 SE) under Q. rotundifolia. At a depth of 30 cm, the contents were lower, ranging from 25.86 g kg⁻¹ (±1.45 SE) to 10.52 g kg⁻¹ (±1.06 SE) respectively.

When the distribution of carbon among the soil particle size fractions was considered, it was observed that the sand fraction concentrated most of the carbon at both depths and in both stands. In the 10 cm layer, values ranged from approximately 26.39 g kg⁻¹ (±1.61 SE) to 19.20 (±0.66 SE) under Q. pyrenaica and Q. rotundifolia. At a depth of 30 cm, the values were approximately 19.53 g kg⁻¹ (±1.97 SE) in Q. pyrenaica and 8.70 g kg⁻¹ (±0.25 SE) in Q. rotundifolia. In the silt fraction, carbon content was substantially lower, ranging in the surface layer from 3.81 g kg⁻¹ (±0.25 SE) in Q. pyrenaica to 2.56 g kg⁻¹ (±0.29 SE) in Q. rotundifolia in the surface layer. At a depth of 30 cm, the values ranged from approximately 2.25 g kg⁻¹ (±0.24 SE) in Q. pyrenaica to 1.39 g kg⁻¹ (±0.25 SE) in Q. rotundifolia. For the clay fraction, carbon content was the lowest of all the analysed fractions. These values ranged from 1.52 g kg⁻¹ (±0.18 SE) to 1.2 g kg⁻¹ (±0.15 SE) in Q. pyrenaica and from 1.54 g kg⁻¹ (±0.07 SE) to 0.61 g kg⁻¹ (±0.05 SE) in Q. rotundifolia.

Soil nitrogen (g kg⁻¹) was significantly influenced by habitat type (Q. rotundifolia vs. Q. pyrenaica) (p < 0.001) and soil depth (p < 0.001) in all fractions analysed (Table 2). Grazing had no significant effect on soil fractionated nitrogen (p > 0.05). None of the interactions between the factors evaluated were significant for this variable (p > 0.05).

The pattern observed for total soil nitrogen (N g kg⁻¹) was similar to that observed for carbon, with differences primarily due to habitat type and soil depth (Figure 4). The means ranged from 2.28 g kg⁻¹ (±0.17 SE) to 0.52 g kg⁻¹ (±0.05 SE) under Q. pyrenaica and Q. rotundifolia, respectively. Among the granulometric fractions, the sand fraction contained the highest concentrations of the N, with averages ranging from of 1.63 g kg⁻¹ (±0.12 SE) under Q. pyrenaica to 0.44 g kg⁻¹ (±0.02 SE) under Q. rotundifolia. In the silt fraction, the values ranged from 0.25 g kg⁻¹ (±0.05 SE) to 0.09 g kg⁻¹ (±0.01 SE), under Q. pyrenaica and Q. rotundifolia, respectively, while in the clay fraction they were lower, ranging from 0.12 g kg⁻¹ (±0.01 SE) in Q. pyrenaica to 0.05 g kg⁻¹ (±0.01 SE) in Q. rotundifolia. In general, N contents decreased with depth in all soil fractions.

Soil carbon stock (Mg ha⁻¹) was significantly influenced by habitat type for total carbon (p < 0.001) and for sand (p < 0.001) and silt (p < 0.05) fractions, while the clay fraction did not show a significant effect (p > 0.05) (Table 3). Grazing and soil depth, alone, had no significant effect on soil carbon stock (p > 0.05). The habitat × soil depth interaction was significant for total carbon stock (p < 0.05) and for sand (p < 0.01) and clay (p < 0.05) fractions, while the other interactions were not significant (p > 0.05).

Soil carbon stock (C, Mg ha⁻¹) varied according to both habitat type and soil depth (Figure 3). In the surface layer, the mean values ranged from 30.42 Mg ha⁻¹ (±2.12 SE) under Q. pyrenaica to 13.25 Mg ha⁻¹ (±2.15 SE) under Q. rotundifolia. Of the granulometric fractions, the sand fraction contained the most carbon, with mean values of 22.51 Mg ha⁻¹ (±2.52 SE) in Q. pyrenaica and 10.61 Mg ha⁻¹ (±2.10 SE) in Q. rotundifolia. In the silt fraction, mean values ranged from 3.54 Mg ha⁻¹ (±0.98 SE) in Q. pyrenaica to 2.42 Mg ha⁻¹ (±0.52 SE) in Q. rotundifolia. The lowest values were found in the clay fraction, ranging from 1.86 Mg ha⁻¹ (±0.32 SE) in Q. pyrenaica to 0.9 Mg ha⁻¹ (±0.28 SE) in Q. rotundifolia. Overall, the highest carbon stocks were observed under Q. pyrenaica, particularly in the sand fraction and in the deeper soil layers.

The results showed that habitat and soil depth were the main determinants of carbon fraction dynamics, whereas the effect of grazing was limited and inconsistent across fractions. Q. pyrenaica exhibited higher C content, particularly in the surface layer (0–10 cm), in terms of both total C and C associated with the sand, silt and clay fractions. This pattern is consistent with a greater potential for carbon accumulation and retention in the Q. pyrenaica stand, although site-specific environmental conditions may also have contributed to the observed differences.

These findings suggests that intrinsic attributes of the habitat, such as greaterlitterfall and litter quality (Zhang et al., 2016; Fonseca and Figueiredo, 2018), root architecture, and intensity of biogeochemical cycling, may play a central role in the formation and stabilisation of carbon in different granulometric fractions (Turrión et al., 2009; Fernández-Getino et al., 2018).

The consistent reduction in content in the 10–30 cm layer, particularly in Q. rotundifolia, highlights the importance of vertical stratification in controlling the distribution of carbon, particularly in terms of total carbon and for the sand and silt fractions. This pattern indicates either lower organic input at depth or greater susceptibility to mineralisation (Jobbágy and Jackson, 2000). Although the absolute values were lower in the clay fraction, the variations associated with the interaction between stand and depth suggest differences in organomineral stabilisation mechanisms between systems (Zhang et al., 2024). These results suggest that the variability in the fractional distribution of carbon is primarily influenced by ecological and structural attributes of Mediterranean oak forests, as well as the vertical organisation of the soil, rather than by pastoral management (Abdalla et al., 2018).

The observed patterns are consistent with the literature on the physical fractionation of soil organic carbon. Studies indicate that most C tends to be associated with macroaggregates, including C enclosed by macroaggregates, while finer fractions (microaggregates and silt + clay) have lower absolute concentrations but greater relative stability (Howlett et al., 2011; Gelaw et al., 2015; Ferreiro-Domínguez et al., 2016). A synthesis of 41 studies by King et al. (2019) showed that macroaggregates account for approximately 83% of the variations in total C stock, with proportionally greater increases in soils with higher carbon contents. The predominance of higher carbon content in the surface fractions associated with macroaggregates, observed in this study, is consistent with this global pattern and reinforces the role of aggregation as a fundamental mechanism for carbon protection.

The absence of a significant effect of grazing on the silt and clay fractions, particularly the latter, which offers the greatest physical-chemical protection, suggests that the management strategy did not encourage sufficient aggregate restructuring to affect the stabilisation of C in the fine fractions. This result is consistent with studies of Mediterranean silvopastoral which show that moderate or rotational grazing does not necessarily result in losses of carbon associated with stable aggregates (Wang et al., 2016; Manetti et al., 2025; Ishaq et al., 2025). Similarly, Lai and Kumar (2020) observed that controlled management systems maintained aggregate stability and carbon associated with macro and microaggregate fractions, even under continuous animal pressure. Therefore, the data suggest that grazing was not a disruptive driver of the structural dynamics of carbon in these areas under a put-and-take grazing system.

Total C content ranged from 10 to 32 g kg⁻¹, which is similar to Turrión et al. (2009) (1.3–18 g kg⁻¹) for different stands of Q. pyrenaica in central-western Spain. This is lower than the values reported by Castro and Fernandez-Nuñez (2014) (6–60 g kg⁻¹) for different Quercus stands in north-eastern Portugal, and substantially lower than the values described by Fonseca et al. (2019) for Q. pyrenaica in north-eastern Portugal (55–100 g kg⁻¹ at 0–10 cm). Similarly, the estimated stocks (13–31 Mg ha⁻¹) fall within the range reported for Mediterranean silvopastoral systems, though they are often lower than the values observed in mature forests or with less disturbed systems. Cardinael et al. (2017) reported values 30–60 Mg ha⁻¹ (0–30 cm) in French agroforestry systems, while Díaz-Pinés et al. (2011) found values of 40–80 Mg ha⁻¹ in Q. pyrenaica stands in central Spain. Lara-Gómez et al. (2020) reported values of 42 Mg ha⁻¹ at depths of up to 30 cm in systems with Q. ilex and Q. suber under traditional management. Conversely, Paustian et al. (2016) described lower values (15–35 Mg ha⁻¹ at 30 cm) were described by Paustian et al. (2016) in temperate pastures with a history of intensive use. Thus, the stocks observed in this study are in an intermediate range, reflecting both soil characteristics and management history.

Regarding nitrogen, all soil layers and fractions had concentrations below 2.3 g kg⁻¹ in ecosystems dominated by Q. pyrenaica, with the exception of total N in the surface layer (0–10 cm). These values are consistent with those reported for Mediterranean Quercus forests. Castro and Fernández-Núñez (2014) reported values of 1.03–8.80 g kg⁻¹ for Q. pyrenaica and 0.61–2.56 g kg⁻¹ for Q. rotundifolia in various locations in Trás-os-Montes (northeastern Portugal), suggesting that the values observed in this study fall within the regional variability previously reported for these ecosystems. In systems dominated by Q. rotundifolia, total N content remained within the range previously reported for Mediterranean Quercus forests, including the values observed for Q. ilex in north-eastern Spain (Canadell and Vilá, 1992) and Q. suber in southern Portugal (Otieno et al., 2011).

In addition to the differences observed between species, the vertical profile of N followed the typical pattern seen in terrestrial ecosystems, with higher concentrations in the surface layer and a progressive decrease down to 30 cm (Wang et al., 2022; Li et al., 2022; Chau et al., 2024). Higher concentrations under Q. pyrenaica (2.2–2.3 g kg⁻¹) than under Q. rotundifolia (1.1–1.2 g kg⁻¹) at the surface suggest a marked effect of stand type on soil N accumulation. This is probably associated with differences in litter quantity and quality, root architecture, and associated microbial activity (Lorenz et al., 2020; Sarıyıldız et al., 2022; Huang et al., 2023; Bashir et al., 2024; He et al., 2023). Overall, these results suggest that the distribution of soil N is influenced by the interaction between habitats (stand-forming tree species) characteristics, edaphic properties, and biogeochemical processes that regulate the decomposition and accumulation of organic matter.

Together, the results support the idea that silvopastoral systems can withstand moderate disturbances fairly well (Mosquera-Losada et al., 2018). The stability of the main carbon compartments in the 0–30 cm layer following grazing. Although the assessment covered a relatively short period, this response may reflect physical protection mechanisms, particularly aggregation, which can help preserve organic matter integrity under disturbance (Sharrow, 2007; Byrnes et al., 2018; Jensen et al., 2019; Meyer et al., 2025).This reinforces, the compatibility between controlled grazing and surface carbon conservation, providing further evidence that well-managed practices can maintain or even sustain the ecosystem service of carbon sequestration (Upson et al., 2016; Teague et al., 2016; Bastani et al., 2023; Guimarães et al., 2023; Lamb et al., 2024). Overall, these results highlights the importance of tree species composition, particularly the dominance of Q. pyrenaica, in controlling soil carbon and nitrogen dynamics in Mediterranean silvopastoral systems. However, given the lack of habitat replication, these differences should be interpreted cautiously, as site-specific factors may also have influenced the observed patterns.

Within the timeframe of this study, soil carbon pools remained relatively stable under the grazing regimes evaluated. However, long-term monitoring is needed to determine whether these systems can consistently sustain their carbon sequestration potential.