Application of Soil Multiparametric Indices to Assess Impacts of Grazing in Mediterranean Forests

Abstract

In this study, the effects of different stocking rates were quantified in three study areas in a Mediterranean forest (Cuenca, Spain) by applying a multiparametric soil quality index (SQI) developed from undisturbed forest soils (>40 years). The main objective was to advance the development and application of multiparametric indices that allow for soil condition assessment. To fulfill this objective, the effectiveness of the developed multiparametric soil quality index (SQI) was analyzed as an indicator of livestock impacts on soil in the Mediterranean forest. The control areas without livestock activity were forest stands of different ages (a thicket forest stand of <30 years; a high-polewood forest stand of 30–60 years; and an old-growth forest stand of >60 years), which were compared with areas subjected to various grazing intensities (areas with permanent livestock passage: a sheepfold that had been inactive for 2–3 years and an active sheepfold; areas with intermittent livestock passage: a bare-soil area, a pine stand and a scrubland). The applied multiparametric soil quality index (SQI) was sensitive to changes in forest ecosystems depending on the stocking rates. However, to obtain greater precision in the assessment of the effects of stocking rates, the multiparametric index was recalibrated to create a new index, the Soil Status Index by Livestock (SSI_L). The correlation between the quality ranges obtained with both indices in different study areas suggests that the SSI_L can be considered a livestock impact reference indicator in Mediterranean forest soils.

1. Introduction

Sheep farming is an important source of wealth in the Cuenca mountain range (Spain), being traditionally practiced alongside other agricultural activities. Most of the region is bounded by extensive areas of livestock practice, mostly associated with areas in which few economic alternatives are found. That is why sheep farming plays a key role in terms of territorial articulation, environment preservation and employment promotion in rural areas in Castilla–La Mancha [1]. Spain represents around 10% of the total sheep and goat livestock in the European Union, and Castilla–La Mancha’s contribution to the national total numbers has risen to 11% according to official agency data provided by the Spanish government.

A large stocking rate (LSR), determined as a measuring unit that compares different livestock categories depending on the type of animal and its age or size (R.D 1053/2022 and R.D 1131/2010), is used to establish the maximum capacity of a farm. Using this approach, the average size of livestock farms is indicated by an LSR of 18.8 per farm in Spain, while in Castilla–La Mancha, it is represented by an LSR of 47.4 per farm (exhibiting wide ranges between autonomous regions). A special mention should be made of Las Majadas (Cuenca district), located in the north of Castilla–La Mancha, where the study areas are situated. This district, with a total pasture area of 3860.19 hectares, has a total LSR of 1126 in extensive livestock practices, which translates to a value of Livestock Loading Unit (LLU = LSR·ha⁻¹) of 0.30 LLU. Of the total value of LLU in the district, 63% corresponds to ovine cattle, 23% to small ruminants, 11.6% to bovine cattle, 1.9% to equine cattle, 0.1% to caprine cattle and 0.4% to wild ungulates such as deer and fallow deer (statistics obtained from the regional forest administration).

Grazing is a biotic factor that affects the ecosystem structure, pasture dynamics [1,2] and microbial soil community, mainly due to soil compaction [3] and changes in the balance between carbon and nutrients [4]. Nevertheless, the response of microbial biomass to grazing is not uniform, sometimes increasing with the intensity of grazing and sometimes not. According to [5], positive effects are commonly visible in high-fertility soil ecosystems, while negative effects are more noticeable in less productive ecosystems.

Even though forest–pasture systems are considered a sustainable alternative for integrated livestock management [6], prolonged grazing could be harmful for soil [7]. Defoliation, trampling and excretion related to intense livestock use of soil [8] are the causes of disturbances affecting its physicochemical properties, such as micro- and macro-nutrient availability (nitrogen, phosphorus and potassium), carbon release from organic matter reservoirs [9] and salinity profile [10], and contribute to soil erosion [11]. In the same way, the composition of soil microorganisms, which take part in recycling nutrients coming from plant wastes and animals, is also altered in forest–pasture systems [5,12].

Multiparametric soil quality indices include physical, chemical and biological soil properties so that the balance can be reflected when certain environmental conditions are met to ensure soil functionality. A valid soil quality index must be sensitive to external disturbances [13] in order to assess the soil state.

The aim of this paper is to advance the application of multiparametric soil indices as indicators of the targeted environmental impacts, specifically impacts caused by livestock practices in Mediterranean forest ecosystems. Particularly, the objectives are (i) to apply a soil quality index (SQI), obtained from Mediterranean forest soils that have been undisturbed for at least the last 20–40 years (five forests located in the Cuenca mountain range, Spain) [14], to soils subjected to varying stocking rate intensities to evaluate sensitivity as an indicator of livestock impacts (active soils, 3-year inactive soils and soils with occasional or null activity), and (ii) to apply a calibration procedure to readjust and obtain a new multiparametric index that confers greater precision when evaluating the impact of grazing activities on forest soils.

2. Materials and Methods

2.1. Study Areas and Experimental Design

The study areas were situated in a forest stand in the Cuenca province (Spain) in a natural park of the Cuenca mountain range named “Ensanche de Las Majadas” (forest code: MUP 133). It is located in East–Central Spain (Figure 1), where the continental Mediterranean climate creates cold winters and mild summers; the mean altitude is around 1440 m, the average rainfall is 934 mm and the mean temperature is 9.10 °C with a range from −3.70 °C to 28.50 °C (Spanish Ministry of Agriculture, Fisheries and Food).

The predominant type of soil is classified as Leptosol [15], with rocky outcrops. Similar conditions can be found across the study areas, except for the bare-soil (BS) and scrubland (Scr) areas, where rocky areas are more frequent.

The control areas’ main forest species is the Spanish black pine (Pinus nigra Arn. ssp. salzmannii), which is sometimes mixed with sporadic individuals of Scots pine (Pinus sylvestris L.) in the pine stand (Ps) areas. The shrub layer´s composition (Scr) is formed of prickly junipers (Juniperus oxycedrus L.), hawthorns (Crataegus monogyna Jacq.), common barberries (Berberis hispanica Boiss & Reut) and common box (Buxus sempervirens L.). Moreover, other species, such as Thymus, Lavandula and Eryngium, appear randomly all over the area. A description of the study areas’ vegetation is shown in

Table 1. Tree and shrub species, vegetation cover density and main characteristics of the study areas in “Ensanche de Las Majadas”.

Area	UTM		Alt. (m)	VEGETATION COVER
	X	Y		TVc (%)	Tree Species	Vc (%)		Shrub Species	Vc (%)		Other Genres		Vc (%)
Ash	585,166	4,459,908	1393	50							Pasture, Thymus, Lavandula, Eryngium		50
Ish	584,979	4,459,104	1446	50									50
BS	585,065	4,458,849	1440	30				Juniperus communis L., Juniperus oxycedrus L., Crataegus monogyna Jacq., Amelanchier ovalis, Berberis hispanica Bois. & Reut., Buxus sempervirens L.	30
Scr	585,822	4,457,380	713	50					30				25
Pst	585,163	4,457,700	1384	70	Pinus nigra Arn. ssp salzmannii, Pinus sylvestris L.	50			50				10
Tfst	586,744	4,458,952	1444	70		50			50				20
Hfst	586,89	4,459,488	1457	70		50			50				20
Ofst	587,651	4,459,046	1444	70		50			50				20

Abbreviations: Ash, active or functional sheepfold; Ish, sheepfold that has been inactive for 3 years; BS, bare soil; Scr, scrubland; Pst, pine stand; Tfst, thicket forest stand; Hfst, high-polewood forest stand; Ofst, old-growth forest stand; UTM, coordinates; Alt., altitude; TVc, total vegetation cover; Vc, vegetation cover. The color intensity corresponds to the type of roof of the zones (three species, green; shrub species, orange; other species, pink).

, along with the dominant species and coverage (Figure 1, Table 1).

Three Spanish black pine stands (Pinus nigra Arn. ssp. salzmannii) of different stand ages and without livestock grazing were selected as the control areas: (i) a thicket forest stand (Tfst) less than 30 years old; (ii) a high-polewood forest stand (Hfst) between 30 and 60 years old; and (iii) an old-growth forest stand (Ofst) over 60 years old. In these control areas, the presence of native wild ungulates is very sporadic, with an LLU of 0.00. At the same time, five areas with livestock activity were selected. Two of these areas are permanently used by livestock: (iv) an active sheepfold (Ash) with an LLU of 45.00, and (v) a sheepfold that has been inactive for at least 2–3 years (Ish) with an LLU of 0.30. Even though the remaining areas exhibit intermittent livestock activity and are habitually used by cattle, they were selected according to their vegetation coverage: (vi) a bare-soil (BS) area; (vii) an area of shrubland–pasture type (Scr); and (viii) a pine stand (Pst). These three areas, located in transhumance routes, have an LLU of 0.30 (

Table 2. Estimation of Livestock Loading Unit (LLU) for the study areas.

Area	Type of livestock	n° Individuals	Cf	LSR	Total LSR	Area (ha)	LLU
Ash	Ovine cattle	1256	0.15	188.40	188.40	4.20	45.00
Ish	Ovine cattle/ Caprine cattle/ Bovine cattle/ Equine cattle	2219	0.15	332.85	1125.90	3860.19	0.30
BS		1	0.15	0.15			0.30
Scr		469	1.00	469.00			0.30
Pst		162	0.40	64.80			0.30
Tfst	Deer/ Fallow deer	17	0.40	6.80	6.80	3964.94	0.00
Hfst							0.00
Ofst							0.00

Abbreviations: Ash, active or functional sheepfold; Ish, sheepfold that has been inactive for 3 years; BS, bare soil; Scr, scrubland; Pst, pine stand; Tfst, thicket forest stand; Hfst, high-polewood forest stand; Ofst, old-growth forest stand; LSR, large stocking rate; Cf, conversion factor. LSR = n° individuals Cf; LLU = LSR·ha⁻¹. Reference: statistics obtained from the regional forest administration, R.D 1053/2022 and R.D 1131/2010.

).

The sampling unit was (10 × 10) m for the sampling plots (total plots = 8), with n = 3 samples (2 × 2) m in each plot, except in the control areas, where n = 2. Each sample was taken from the first 15 cm of soil, removing leaf litter, and was composed of 6 subsamples in order to minimize spatial variability [16,17]. Sampling took place for a whole annual cycle (spring, summer, autumn and winter) starting in autumn 2016. A total of 84 samples were taken, with 6 samples from the intense livestock activity areas (3 samples × 2 plots), 9 samples from the intermittent livestock activity areas (3 samples × 3 plots) and 6 samples from the control areas (2 samples × 3 plots). Across the 4 seasons, this constituted a total of n = 168 samples (two replicas included).

2.2. Parameters Analyzed in This Study

A series of physicochemical, microbiological and enzymatic soil parameters were analyzed: gravimetric moisture (M, %); pH using a pH meter; electric conductivity (EC, mS·m⁻¹) (Navi Horiba model); total organic carbon (TOC, %) [18]; nitrogen (N,%) [19]; phosphorus (P, ppm) [20]; basal respiration (BR, μgC-CO₂·g⁻¹·day⁻¹) [21]; microbial biomass carbon (MBC, μgC·g⁻¹) using the fumigation–extraction method proposed by Vance et al., 1987, and adapted in [22]; dehydrogenase enzyme activity (DHA, μmol (INTF) ·g⁻¹·h⁻¹) [22]; alkaline phosphatase (APA, μmol (PNF)·g⁻¹·h⁻¹) and β-glucosidase (β-GLU, μmol (PNF)·g⁻¹·h⁻¹) as determined following the Tabatabai and Bremmer methods [23]; and urease enzymatic activity (UA, μmol (N-NH₄⁺)·g⁻¹·h⁻¹) [24].

2.3. Soil Quality Index (SQI)

A multiparametric soil quality index (SQI) was applied, which was developed based on unaltered Mediterranean forest soils that had been unmanaged for at least for 20 to 40 years [14]. These soils were mainly obtained from 5 forests located in the Cuenca mountain range (Spain), including the study area “Ensanche de Las Majadas”, where the main species is the Spanish black pine (Pinus nigra Arn. ssp. salzmannii). The index is shown in Equation (1), which was obtained from a statistical analysis of 12 physicochemical, microbiological and enzymatic variables. After subsequently using the method of principal component analysis (PCA), the selected representative variables were moisture (M), pH, total organic carbon (TOC), microbial biomass carbon (MBC), and alkaline phosphatase (APA) and β-glucosidase (βGLU) enzymatic activities.

$$\begin{gathered} \mathrm{S}\mathrm{Q}\mathrm{I}=0.576·\left[0.489·\left(\frac{1}{1+{\left(\frac{1308 - 16.31}{\mathrm{M}\mathrm{B}\mathrm{C} - 16.31}\right)}^2}\right)+0.459·\left({\mathrm{e}}^{ - \left(\frac{(\mathrm{M} - {39)}^2}{2·{11.2}^2}\right)}\right)+0.445·\left(\frac{1}{1+{\left(\frac{50 - 19.6}{\mathrm{T}\mathrm{O}\mathrm{C} - 19.6}\right)}^2}\right)\right]+0.228 \\ ·\left[ - 0.602·\left(\frac{1}{1+{\left(\frac{121.9}{\mathrm{A}\mathrm{P}\mathrm{A}}\right)}^{1.7}}\right)+0.510·\left(\frac{1}{1+{\left(\frac{197.2}{\mathsf{\beta }\mathrm{G}\mathrm{l}\mathrm{u}}\right)}^{1.7}}\right)\right]+0.196·\left[0.831·\left({\mathrm{e}}^{ - \left(\frac{(\mathrm{p}\mathrm{H} - {6)}^2}{2·{0.59}^2}\right)}\right)\right] \end{gathered}$$

(1)

The values for soil quality in each study area were obtained by applying Equation (1) and standardized to between 0 and 1 to establish environmental quality ranges. Every result was linked to the maximum and minimum values of the whole group of data [25]. This allowed for an accurate comparison between different situations or activities in the environment.

A quartile division was carried out to establish different soil quality levels to be taken as the reference: values of 0.00–0.25 indicate low quality, 0.26–0.50 indicate low–medium quality, 0.51–0.75 indicate medium–high quality and 0.76–1.00 indicate high quality (

Table 3. Soil quality value ranges.

Range	Quality
0.76–1.00	A	High
0.51–0.75	B	Medium–high
0.26–0.50	C	Medium–low
0.00–0.25	D	Low

).

2.4. Soil Multiparametric Index Development (SSI_L)

A similar methodological procedure to the one used in development of the SQI was applied to the study areas (control areas and areas with soils of various stocking rate intensities) to calibrate a new index that could evaluate livestock impacts with greater precision and, thus, better reflect soil condition due to livestock activity; this would allow for an analysis of which activities are sustainable and which are not. This new index was defined as the Soil Status Index by Livestock (SSI_L), following the normalization and selection process described in [14], which consists of (i) selecting representative parameters by means of consecutive principal component analysis (ACP); (ii) transforming and normalizing data using standardized functions for each parameter according to its contribution to environmental quality; and (iii) analyzing and combining values into a model, while bearing in mind each component’s value and the selected parameter weight. Finally, a mathematic function was obtained, which was a combination of the selected parameters into a multiparametric soil condition index.

The steps were as follows: (I) A group of 12 physicochemical, microbiological and enzymatic soil variables were analyzed: M, pH, EC, TOC, N, P, BR, MBC, DHA, UA, APA and β-GLU. (II) A first principal component analysis, 1ACP, was carried out with all the variables, and those with higher eigenvalues (>1.00) were selected as the principal components (PC). For each of the selected principal components, those with the maximum eigenvalue and those with an eigenvalue within 90% of the maximum were chosen as the representative parameters [26]. (III) To discard correlated variables, a second principal component analysis, 2ACP, was carried out using only those parameters that were selected beforehand. (IV) A third and last principal component analysis, 3ACP, was carried out to obtain the coefficients related to the selected parameters included in the index. These variables were standardized by using the function “more is better” [27] and a Gaussian-type function [28]. (V) Lastly, the index was defined by a linear combination of those transformed and weighted values and given as a result of the SSI_L equation (Soil Status Index by Livestock).

2.5. Statistical Analysis

An analysis of variance was applied using general linear models (GLMs) to (i) characterize the physicochemical, microbiological and enzymatic parameters between different study areas and (ii) study the soil quality index (SQI)’s sensitivity in order to analyze the effects of stocking rate intensity, season and the interactions between these two on the SQI results.

A consecutive principal component multifactorial analysis was carried out to calibrate the new Soil Status Index by Livestock (SSI_L), and its response significance to livestock use was assessed with a GLM. Regression models were applied to study the relations between the SQI and SSI_L.

Fisher’s Least Significant Difference (LSD) method (95% confidence interval) was used, with p < 0.05, and applied to cases where a significant F value was obtained. The software used for the statistical analysis was Statgraphics centurion XVI.

3. Results

3.1. Physicochemical, Microbiological and Biochemical Characterization of Soils

Table 4. Mean values and standard deviations for each parameter in each study area (n = 168).

Parameters I	Physicochemical						Microbiological		Enzymatic Activities
	TOC	N	P	M	pH	EC	BR	MBC	DHA	UA	APA	β-GLU
	(%)	(%)	(ppm)	(%)		(mS · m−1)	(μg C-CO2·	(μg C · g−1)	(μmol(INTF)	(μmol(N-NH4+)	(μmol(PNF)	(μmol(PNF)
Area II							g−1· day−1)		· g−1 · h−1)	· g−1 · h−1)	· g−1 · h−1)	· g−1 · h−1)
Ash	5.61 ± 0.23 b	1.15 ± 0.06 a	191.78 ± 15.84 a	32.33 ± 1.26 b	7.55 ± 0.08 cd	31.66 ± 1.67 a	126.55 ± 9.93 ab	2142.25 ± 120.44 a	0.03 ± 0.00 cd	27.65 ± 1.48 a	116.39 ± 7.00 a	34.77 ± 1.90 d
Ish	5.07 ± 0.23 b	1.01 ± 0.06 abc	187.10 ± 15.84 a	23.95 ± 1.26 d	7.72 ± 0.08 abc	20.45 ± 1.67 bc	48.40 ± 9.93 ef	1073.11 ± 120.44 bc	0.03 ± 0.00 bc	21.78 ± 1.48 b	94.36 ± 7.00 cd	28.28 ± 1.90 e
BS	5.26 ± 0.23 b	0.88 ± 0.06 c	48.66 ± 15.84 b	19.67 ± 1.26 e	7.92 ± 0.08 a	14.55 ± 1.67 d	35.41 ± 9.93 f	821.48 ± 120.44 c	0.04 ± 0.00 a	19.22 ± 1.48 b	84.80 ± 7.00 de	37.15 ± 1.90 cd
Scr	5.18 ± 0.23 b	0.87 ± 0.06 c	37.98 ± 15.84 b	24.25 ± 1.26 d	7.91 ± 0.08 ab	16.05 ± 1.67 cd	31.43 ± 9.93 f	792.68 ± 120.44 c	0.04 ± 0.00 ab	20.53 ± 1.48 b	112.75 ± 7.00 abc	32.22 ± 1.90 de
Pst	6.77 ± 0.23 a	1.09 ± 0.06 ab	54.11 ± 15.84 b	38.26 ± 1.26 a	7.33 ± 0.08 d	31.05 ± 1.67 a	99.15 ± 9.93 bc	1041.76 ± 120.44 bc	0.03 ± 0.00 de	13.89 ± 1.48 c	115.08 ± 7.00 ab	51.73 ± 1.90 b
Tfst	5.29 ± 0.28 b	0.86 ± 0.07 c	43.44 ± 19.39 b	27.65 ± 1.55 cd	7.66 ± 0.10 bc	24.92 ± 2.05 b	67.23 ± 12.16 de	844.04 ± 147.51 bc	0.03 ± 0.00 e	12.34 ± 1.81 c	67.53 ± 8.57 e	42.70 ± 2.33 c
Hfst	5.19 ± 0.28 b	0.93 ± 0.07 bc	43.18 ± 19.39 b	30.51 ± 1.55 bc	6.74 ± 0.10 e	17.67 ± 2.05 cd	88.86 ± 12.16 cd	1028.94 ± 147.51 bc	0.01 ± 0.00 f	12.16 ± 1.81 c	77.62 ± 8.57 de	57.45 ± 2.33 ab
Ofst	6.71 ± 0.28 a	1.08 ± 0.07 ab	46.31 ± 19.39 b	37.40 ± 1.55 a	5.97 ± 0.10 f	17.24 ± 2.05 cd	142.21 ± 2.16 a	1249.34 ± 147.51 b	0.01 ± 0.00 f	12.12 ± 1.81 c	93.21 ± 8.57 bcd	59.86 ± 2.33 ab

^I TOC, total organic carbon; N, total nitrogen; P, phosphorus; M, moisture; pH, soil acidity; EC, electrical conductivity; BR, basal soil respiration; MBC, microbial biomass carbon; DHA, dehydrogenase activity; UA, urease activity; APA, phosphatase activity; β-GLU, β-glucosidase activity. ^II. Ash, active or functional sheepfold; Ish, sheepfold that has been inactive for 3 years; BS, bare soil; Scr, scrubland; Pst, pine stand; Tfst, thicket forest stand; Hfst, high-polewood forest stand; Ofst, old-growth forest stand. ^{a, b, c, d, e, f} Homogeneous subgroupings according to the results of the analysis of variance.

shows all 12 edaphic variables analyzed (TOC, N, P, M, pH, EC, BR, MBC, DHA, UA, APA and β-GLU) for the study areas. It is observed that physicochemical parameters, such as TOC, form two uniform subgroups. Higher values (more than 6.70) are registered in Pst and Osh, while the rest of the study areas show lower values (5.61 measured in Ash). On the contrary, the highest value of N appears in Ash (1.15) without any significant differences between Pst, Ofst and Ish (1.09, 1.08 and 1.01, respectively). The remaining areas present significantly lower values, with Scr (0.87) and BS (0.88) being those areas with the lowest N concentrations. Ash and Ish show significantly higher p values than the other areas (191.78 and 187.10, respectively), while Scr is registered as having the lowest value (37.98). As for moisture, the areas with the highest values are Pst (38.26) and Ofst (37.40), while BS is the area with the lowest value (19.67). Lower pH values are found in Ofst (5.97), followed by Hfst (6.74), while BS (7.92), Scr (7.91) and Ish (7.72) are the areas with the highest pH values. The EC values are higher in Ash (31.66) and Pst (20.45), while the values in the other areas range between these results, with the lowest value found in BS (14.55). Both BS and Scr show the lowest values of MBC and BR, while Ofst and Ash have the highest values (Table 4). In relation to enzymatic activities, Ash shows significantly higher values in terms of UA (27.65) and APA (116.39), while the control areas have the lowest values for these two parameters. The lowest value of DHA is also found in the Hfst and Ofst control areas (0.01), while BS and Scr have a significantly higher value (0.04). The remaining areas are characterized by an intermediate value (0.03). However, β-GLU shows the highest values in the Ofst (59.86) and Hfst (57.45) control areas and intermediate values in Pst (51.73) and Tsft (42.70), while the rest range between these values and the value found in Ish (28.28), which is at the lowest end for that parameter.

3.2. Sensitivity Analysis of the Applied Soil Quality Index (SQI)

The results in

Table 5. Significance levels of the factors intensity zone (IZ), season (S) and their interaction (IZ × S) and their effects on the SQI value.

	GLM
	Fp			Model
Variable	IZ	S	IZ × S	R2	Fp
SQI	20.15 ***	96.91 ***	3.90 ***	79.20	16.70 ***

Model fit level [F: Snedecor’s F-distribution; R²: coefficient of determination; p < 0.001 (***); n = 168]. General lineal model, GLM.

show that the stocking rate, the season and the interactions between these two variables have a significant influence (p < 0.001) on the SQI values, explaining 79.20% of the total variability.

The mean SQI values are the highest in the Ofst (0.42) control area, while the lowest values are registered for Ish, Scr and BS with intermittent livestock passage (0.24, 0.19 and 0.19, respectively). The other study areas, including the active sheepfold area (0.32), show intermediate values (Figure 2,

Table 6. Mean values of the soil quality index (SQI) in each study area.

Area	SQI
Ash	0.32	±	0.02	bc
Ish	0.24	±	0.02	de
BS	0.19	±	0.02	f
Scr	0.19	±	0.02	ef
Pst	0.25	±	0.02	d
Tfst	0.27	±	0.02	cd
Hfst	0.35	±	0.02	b
Ofst	0.42	±	0.02	a

Ash, active sheepfold; Ish, inactive sheepfold; BS, bare soil; Scr, scrubland; Pst, pine stand; Tfst, thicket forest stand; Hfst, high-polewood forest stand; and Ofst, old-growth forest stand (n = 168, units shown on the planes of the axes). ^{a, b, c, d, e, f} Homogeneous subgroups, p < 0.05.

).

3.3. New Soil Multiparametric Index Development (SSI_L)

The first principal component analysis was carried out (1PCA) including the full set of variables considered (12); the results show that three of the components have an eigenvalue of ≥1, explaining 66.36% of the cumulative variance (CV) in the data (

Table 7. Results of a principal component analysis (PCA) performed with the full set of parameters: eigenvalues from the first three principal components, percentage variance explained (EV), cumulative variance percentage (CV) and corrected explained variance (CEV).

Principal Component Analysis	1PCA			2PCA			3PCA
Principal Component	1PC1	1PC2	1PC3	2PC1	2PC2	2PC3	3PC1	3PC2
Eigenvalue	4.46	1.94	1.55	3.25	1.69	1.03	1.55	1.11
EV	37.18	16.21	12.96	40.66	21.08	12.91	51.61	37.11
CV	37.18	53.39	66.36	40.66	61.74	74.65	51.61	88.72
CEV							58.17	41.82
Parameters a
TOC	0.36	−0.09	0.13	0.40
N	−0.08	−0.33	0.48
P	0.07	0.55	0.14
M	0.36	−0.24	0.20	0.43			0.72
pH	−0.35	0.09	0.42	−0.42
EC	0.31	0.29	0.00
BR	0.35	0.05	0.09	0.43
MBC	0.38	0.20	0.20	0.43			0.69
DHA	−0.23	0.01	0.54			0.73
UA	−0.04	0.57	0.04		0.67			0.92
APA	0.33	−0.03	0.37
B-GLU	0.29	−0.23	−0.19

^a Bold values correspond to larger eigenvectors (>90% of the maximum weight per PC). TOC, total organic carbon; N, total nitrogen; P, phosphorus; M, moisture; pH, soil acidity; EC, electrical conductivity; BR, basal soil respiration; MBC, microbial biomass carbon; DHA, dehydrogenase activity; UA, urease activity; APA, phosphatase activity; β-GLU, β-glucosidase activity.

, Figure 3a). The selected parameters for each component are those with a higher weight (in bold and an intense shade) and that are within the 90% range (in bold and a less intense shade). The scatter plot (Figure 3b) shows the distribution areas according to the parameters.

When a second principal component analysis (2PCA) was carried out, a strong correlation between the highest-weighted variables was revealed (Figure 3c); thus, it was decided to select those variables because they are included in the soil quality index (SQI) used as the reference.

In the third principal component analysis (3PCA), the cumulative variance was corrected in 100% of the cases (CEV). The corrected explained variance (CEV) together with each parameter´s weight is listed in the developed model for the new index, the SSI_L (Equation (2)).

The selected parameters for each consecutive ACP are shown in Figure 3d.

Moisture (M), microbial biomass carbon (MBC) and urease activity (UA) were the parameters selected to be part of the SSI_L. In

Table 8. Values of the constants of each standardization equation and correlation coefficient for each of the parameters that make up the SSI_L.

Variables a or Constants b	Op	σ	B	L	m	Standarization Equation	r
M (%)	28.88	13.79				y = e − ((M − 28.88)2)/(2 × 13.792))	0.99
MBC (μg C g−1)			1136.12	0.01	1.8	y = 1/(1 + (1136.12) / (MBC))1.8)	0.99
UA (μmol N-NH4+ g−1 h−1)			18.21	1.33	1.8	y = 1/(1+ (18.21 / UA)1.8)	0.98

^a M, moisture; MBC, microbial biomass carbon; UA, urease activity; ^b Op, optimal value; σ, standard deviation; B, critical value; L, lower value; m, slope of the equation.

, the normalization equation for each variable is displayed, together with the adjustment factor, critical value, optimal value, standard deviation and smallest value.

The SSI_L is the result of the weighted summation of the normalized values of moisture (M), microbial biomass carbon (MBC) and urease activity (UA), following the model shown in Equation (2).

$$\begin{array}{l}{\mathbf{S}\mathbf{S}\mathbf{I}}_{\mathbf{L}}=0.58·[0.72·\left({\mathrm{e}}^{ - \left(\frac{{\left(\mathrm{M} - 28.88\right)}^2}{(2 \times {13.79}^2)}\right)}\right)+0.69\\ ·\left(\frac{1}{1 + {\left(\frac{1136.12}{\mathrm{M}\mathrm{B}\mathrm{C}}\right)}^{1.8}}\right)]+0.42\\ ·\left[0.92·\left(\frac{1}{1 + {\left(\frac{18.21}{\mathrm{U}\mathrm{A}}\right)}^{1.8}}\right)\right]\end{array}$$

(2)

3.4. Sensitivity Analysis of the New Multiparametric Index, SSI_L

The results obtained from fitting the linear statistical models to relate the SSI_L to the selected factors (IZ, intensity zone; S, season; and IZ × S, interaction) show the high sensitivity of the SSI_L, with a confidence level > 95.0%, with significant differences between areas with different livestock grazing intensities under a seasonal influence (

Table 9. Significance levels of the factors intensity zone (IZ) and season (S) and their interaction (IZ × S) and their effects on the SSI_L value.

	GLM
	Fp			Model
Variable	IZ	S	IZ × S	R2	Fp
SSIL	11.38 ***	38.70 ***	5.81 ***	70.64	10.55 ***

Model fit level [F: Snedecor’s F-distribution, R²: coefficient of determination; p < 0.001 (***); n = 168].

).

The SSI_L shows the highest value in the area with the highest load intensity, Ash (0.75), and the lowest values in intermittent cattle passages, Pst (0.49), as well as in the control zones, Hfst (0.59), Ofst (0.56) and Tfst (0.56), leaving an intermediate group with no significant differences (Figure 4,

Table 10. Mean values and standard deviations of the Soil Status Index by Livestock (SSI_L) in the study areas.

Area	SSIL
Ash	0.75	±	0.02	a
Ish	0.65	±	0.02	b
BS	0.62	±	0.02	bc
Scr	0.58	±	0.02	c
Pst	0.49	±	0.02	d
Tfst	0.56	±	0.03	cd
Hfst	0.59	±	0.03	bc
Ofst	0.56	±	0.03	cd

Ash, active sheepfold; Ish, inactive sheepfold; BS, bare soil; Scr, scrubland; Pst, pine stand; Tfst, thicket forest stand; Hfst, high-polewood forest stand; and Ofst, old-growth forest stand (n = 168, units shown in the plane of the axes). ^{a, b, c, d} Homogeneous subgroups, p < 0.05.

).

3.5. Correlation between SQI and SSI_L

The best correlation between the two indices was obtained with a simple linear regression model (p < 0.001), which explains 84.31% of the variability (

Table 11. Adjustment equation and correlation between SQI and SSI_L.

	Model
	Fp	R2 (%)	SEE
SQI = 0.4531 · SSIL	897.69 ***	84.31	0.12

Model fit level, SQI variable and adjustment level of the model [F: Snedecor’s F-distribution, R²: coefficient of determination, SEE: standard error of the estimate; all models are significant; p < 0.001 (***); n = 168].

). The correlation between the two indices is positive and relatively strong (correlation coefficient of 0.92), although soil quality is reduced by 45% (as quantified by the SQI) when the effect of livestock activity on soil increases (as quantified by the SSI_L).

The obtained values for soil quality (SQI) and soil condition by livestock use (SSI_L), which were standardized between 0 and 1, allowed a comparison between different situations or actions in the environment according to the different soil quality ranges taken as the reference, as shown in Table 3. In this way, a zone classification was established based on soil quality and livestock usage activity (

Table 12. Classification of soil areas based on quality or livestock activity levels established using the Rank Soil Quality Index (RSQI) and the Rank Soil Status Index by Livestock (RSSI_L).

Area	RSQI	Quality	RSSIL	Status by livestock
Ash	0.49 ± 0.03 bc	medium–high	0.64 ± 0.03 a	medium–high
Ish	0.36 ± 0.03 de	medium–low	0.52 ± 0.03 b	medium–high
BS	0.28 ± 0.03 f	medium–low	0.48 ± 0.03 bc	medium–low
Scr	0.29 ± 0.03 ef	medium–low	0.43 ± 0.03 c	medium–low
Pst	0.38 ± 0.03 d	medium–high	0.32 ± 0.03 d	medium–low
Tfst	0.42 ± 0.03 cd	medium–high	0.41 ± 0.03 cd	medium–low
Hfst	0.55 ± 0.03 b	medium–high	0.44 ± 0.03 bc	medium–low
Ofst	0.67 ± 0.03 a	high	0.40 ± 0.03 cd	medium–low

Ash, active sheepfold; Ish, inactive sheepfold; BS, bare soil; Scr, scrubland; Pst, pine stand; Tfst, thicket forest stand; Hfst, high-polewood forest stand; and Ofst, old-growth forest stand. ^{a, b, c, d} Homogeneous subgroups, p < 0.05.

, Figure 5).

The control areas (Ofst, Hfst and Tfst) and the Pst and Ash areas show medium–high quality according to the SQI range. The other areas are classified as medium–low quality. However, when taking into account the Soil Status Index by Livestock (SSI_L), it can be seen that the Ash and Ish areas fall into the medium–high range, which is different from the other areas and is possibly attributable to these areas’ greater livestock activity.

4. Discussion

According to the physicochemical parameters chosen for this study, micro- and macro-nutrients vary significantly between the study areas. Carbon (TOC) is higher in the control areas and mature stands (Ofst and Pst, respectively), where vegetation coverage contributes to the organic carbon accumulation in the superficial top layers of the soil (0.00–0.05 m) [29]. In addition, in the pine stand (Pst), light grazing enhances primary production and the nutrient cycle in such a pasture ecosystem, as mentioned by [30]. A similar observation was made for N, which also increases in areas with higher livestock intensity (Ash and Ish), where P shows the highest results as well. Livestock depositions and soil compaction due to trampling in high-intensity livestock areas elevate N and P concentrations, although C availability depends more on the stand age and crown type [9,31,32,33]. However, this study did not take into account the compaction or alteration of the bulk density of the soil, as it was not considered a determining factor in forest areas with extensive livestock. In addition, similar studies on soil responses to animal traffic show inconclusive results since their effect depends on various factors, such as the type of livestock, soil moisture and the type of exploitation (intensive or extensive) [34,35]

The relationship between microbiological parameters and enzymatic activities shows that the abundance of microbial biomass carbon (MBC) might be due to the accumulation of depositions coming from stocking rate in the active sheepfold (Ash); stockpiled organic wastes like litterfall in mature stands (Ofs); or both circumstances, as happened in the pine stand areas (Pst). The high concentration of soil microorganisms detected might be related to organic matter accumulation coming from animal excretion or organic wastes in a specific area. This phenomenon might be considered an incentive for microorganism activity and development, playing a key role in organic matter decomposition and improvement in soil health in various ecosystems of the world according to the study carried out in [36].

The high sensitivity shown by enzymatic activities when disturbances happen makes them excellent alternative indicators of soil dynamics [37]. Accordingly, the high values obtained for urease (UA) and phosphatase (APA) activities in Ash and Ish, which are associated with the nitrogen (N) and phosphorus (P) cycles, respectively, indicate important metabolic activity related to these biogeochemical cycles, in addition to the associated high stocking rates [38]. The SQI’s sensitivity was tested to determine its usefulness in evaluating the different areas and their livestock activity. The control areas, where livestock activity is not carried out, show a superior value of soil quality. This is believed to be caused by the essential properties of these unaltered areas in terms of vegetation cover, such as in the case of the old-growth forest stand area (Ofst) with a high contribution of litterfall. Such vegetation remains from pine coverage perform an important role in providing organic matter, preventing soil erosion and increasing the moisture content, which is reflected in the improvement in the mineralization process, which leads to higher nutrient liberation [17].

Moreover, the SQI is proven to be a useful indicator of activity or disturbance in soil functionality. In the intermittent grazing areas, despite the livestock activity, positive effects are observed in high-fertility ecosystems (as in the standard pine forest with intermediate livestock passage, Pst), while negative effects are commonly seen in less productive ecosystems (bare soil (BS), scrubland (Scr) and inactive sheepfold (Ish)) when stocking rates are medium [5].

Therefore, the SQI not only allows for the assessment of soil quality, but also helps understand how livestock activity influences ecosystems with different fertility levels. The results presented in this paper support other studies showing that grazing activities have a positive influence; however, this influence fluctuates in accordance with the involved ecosystems’ characteristics, which suggests that the sensitivity of the SQI allows it to distinguish different impacts depending on the ecosystem conditions [39]. Nevertheless, in the active sheepfold (Ash), a less productive area with a high stocking rate, a high value of the SQI is also obtained. These results might be caused by the significant contribution of organic matter from animal depositions, consequently triggering enzyme and microorganism activation [40,41].

Despite the SQI’s validation, when applied to assess the quality of unaltered soil with mature vegetation coverage, the contradictory data obtained in areas with a high-intensity stocking rate lead to the development of a new index.

The validity and sensitivity of the new Soil Status Index by Livestock (SSI_L) were evaluated using the procedure established in [14], with altered and unaltered soils, to objectively quantify soil conditions against disturbances such as intense livestock activity. Significant differences were observed between the areas with high and intermediate stocking rates and the control areas, suggesting that the SSI_L is capable of effectively quantifying disturbances caused by livestock activity. The highest value of the SSI_L was observed in the active sheepfold, followed by the areas with an intermediate stocking rate, while the lowest values correspond to the control areas [16].

The positive correlation between the normalized quality ranges obtained with both indices (SQI and SSI_L) reinforces the usefulness of the SQI in assessing soil quality in natural ecosystems and the SSI_L in assessing quality disturbances caused by varying intensities of stocking rates. In addition, when comparing the SSI_L with the SQI, different patterns according to the established ranges are revealed, providing an exhaustive comprehension of soil quality’s and livestock activity´s impacts on the environment [42,43].

The applied methodology seems to be sensitive and robust in assessing livestock activity´s impacts on soil quality. These findings highlight the importance of understanding the relationships between livestock activity and soil quality and how using indices such as the SQI and SSI_L can provide valuable information for soil management and livestock activity decision-making processes [44].

5. Conclusions

The SQI (soil quality index) appears to be a versatile tool, as it is sensitive to changes in soil functionality even in the face of impacts such as grazing. Its applicability to undisturbed forest ecosystems and its ability to detect disturbances and differences in productivity in various ecosystems make it a valuable tool.

The SSI_L (Soil Status Index by Livestock) is presented as a specific and objective indicator to quantify the disturbance caused by the stocking rate. Its ability to objectively measure the impact of livestock activity on soil provides valuable information on soil sustainability.

The significant correlation between the two indices reinforces their usefulness for measuring soil sustainability in the face of disturbances, suggesting that the use of both indices can provide a more complete and accurate assessment of soil health in grazed forest ecosystems, and that they are useful tools for decision-making in forest management. However, more studies would be needed to look for new interactions between other less studied soil parameters, such as density and resistance to penetration, that could be tailings in intensively grazed ecosystems.

The results obtained in this study encourage us to continue advancing the objective quantification of environmental impacts through the application of multiparameter indices as indicators, which can facilitate informed decision-making on sustainable forestry practices and soil restoration.