Plant-soil interactions in response to grazing intensity in a semi-arid ecosystem from NE Spain

Abstract Livestock grazing is an important element in ecosystem regulation since it may affect essential ecosystem functions, such as nutrient acquisition, organic matter decomposition, or litter accumulation in the soil. Overgrazing can threaten the conservation of ecosystems through excessive defoliation of plants or trampling. On the contrary, moderate grazing can have benefits on ecosystem dynamics by favoring nutrient cycling or the soil microbial activity. The aim of this study was to analyze these effects in a semi-arid Mediterranean shrubland located in NE Spain. We established six study sites including three grazing intensities, where we sampled vegetation biomass and soil properties: nitrogen content, microbial biomass, water infiltration capacity, porosity, and gypsum content. These parameters were included in a plant-soil interaction model tested through Structural Equation Modeling. Grazing had a direct negative effect on plant biomass (p < 0.01) and water infiltration capacity (p < 0.05) affecting soil nitrogen content (p < 0.001) and microbial biomass (p < 0.5), respectively. Infiltration capacity and porosity were primary drivers of plant biomass (p < 0.05, both cases), and plant biomass was the main contributor to the soil nitrogen pool. Microbial biomass was dependent on infiltration capacity (p < 0.05), porosity (p < 0.01), and nitrogen (p < 0.01). Grazing directly or indirectly affected the functioning of the ecosystem through effects on plant and soil attributes, which may result in changes in plant growth, litter decomposition, or plant nutrient acquisition. This study revealed that moderate grazing can maintain optimal ecosystem features and prevent ecosystem degradation.


Introduction
Grazing and demand for its by-products are expected to increase by the middle of the 21st century (World Bank 2008), which is a potential driver of desertification in resource-poor ecosystems, such as drylands (FAO 2019), which occupy 40% of Earth's surface and harbor 40% of the human population (Reynolds et al. 2007). Among their ecosystem services, drylands supply fiber and food, contribute to the global economy with 30% of cultivated plants and livestock breeds adapted to arid conditions (Lee, Schaaf and UNESCO 2008), provide forage production and improve the preservation and resource status of soil, while maintaining biodiversity (Gratani et al. 2013).
Semi-arid shrublands are a type of dryland widely spread throughout the Mediterranean region, characterized by singular and diverse plant communities generally distributed in fertile patches of shrubs alternating with bare soil interspaces poor in nutrients and microorganisms (Pugnaire et al. 1996). Fertile patches act as ecosystem engineers providing important services, such as facilitation of seed germination or seedling establishment (Pueyo et al. 2016;Foronda et al. 2019), promotion of microbial growth (Garc ıa et al. 2002;Goberna et al. 2007;Maestre et al. 2009), prevention of erosion (Casermeiro et al. 2004), nutrient enrichment of soil (Maestre et al. 2009), or improvement of the water and light regime in soil (De Dato et al. 2009).
Overall, grazers impact plant biomass by browsing (Tessema et al. 2011) and changing the floristic composition, richness, and arrangement of species in the plant community (Angassa 2014). They affect soil nutrient status by modifying the C:N ratio of plant tissue (Baron et al. 2001) and reducing litter fall and thus nutrient returns to soil (Sch€ onbach et al. 2011), which may affect the biomass and activity of microorganisms (Prieto et al. 2011). Grazers may affect the hydro-physical properties of soil by trampling (Wang and Wesche 2016), which impacts the form and stability of soil aggregates and alters the pore size and infiltration capacity of the soil (Taboada, Rubio, and Chaneton 2011;Moret-Fern andez et al. 2021). Grazers also redistribute nutrients by excreta and urine, which are rich in easily decomposable N (Rossignol, Bonis, and Bouzill e 2006), and by movements at the landscape scale (Bailey andBrown 2011, Tessema et al. 2011).
Moderate grazing has been shown to act as a positive driver of ecosystem processes in these arid ecosystems, promoting important ecosystem functions, such as plant growth, nitrogen, and organic carbon stocks, microbial activity or soil hydro-physical properties (Tessema et al. 2011;Liu et al. 2012;Moret-Fern andez et al. 2021). In contrast, overgrazing can have substantial negative effects on ecosystem functions (Cardinale et al. 2012), making the ecosystems more vulnerable to global change and desertification.
Whereas most studies have traced particular and direct effects of livestock grazing on ecosystem properties, such as nutrient status or microbial growth (Baron et al. 2001;Prieto et al. 2011), there is a lack of understanding of a general mechanism explaining the direct and indirect interactions between livestock grazing, vegetation and soil functions in Mediterranean shrublands, and what the relative importance of these interactions is in the whole system. The aim of this study was to elucidate this understanding by relating livestock grazing to plant biomass and soil biological and hydro-physical properties, i.e., resource status, microbial biomass, porosity, and infiltration capacity. In this study, we hypothesized that livestock act on these shrublands through two major pathways: (1) changes in the plant biomass by plant tissues removal and (2) changes in hydro-physical properties of soil by trampling and that (3) these effects affect the nitrogen and water accessibility for plants and microorganisms, with repercussions on their growth, and thus their biomass. Specifically, we hypothesized that (1) grazing has a direct negative effect on plant biomass, linked to the browsing animals (Harrison and Bardgett 2004), which (2) indirectly reduces the nitrogen content in the soil since less litter enters the soil. We also wanted to prove that (3) grazing has a direct negative effect on the soil hydro-physical properties (i.e., porosity and infiltration capacity) potentially caused by the trampling of animals (Moret-Fern andez et al. 2021), that (4) would indirectly reduce plant and microbial biomass through changes in soil structure and water availability (Wen, He, and Zhang 2016).

Study area
The study area was located in the Middle Ebro Valley, NE Spain (41 26 0 0 00 , 0 44 0 0 00 W), where the climate is the semi-arid Mediterranean. The annual mean temperature is around 14.9 C and annual precipitation ranges from 300 to 400 mm À1 (SIAR, Ministerio de Agricultura, Pesca, Alimentaci on y Medio Ambiente). Altitude ranges from 200 to 600 m a.s.l. Soils are gypsiferous, with low content of organic matter (1.94%) and basic pH (7.84), with a sandy loam texture (USDA classification), and classified as Typic Haplogypsids and Gypsic Haplosalids (Moret-Fern andez et al. 2021).

Methodological procedure
Field surveys and soil sampling We set up six independent study sites, each one including three different grazing intensities (low, moderate, and high) on 14 April 2017. At each level of grazing intensity, we marked out five plots of 1 m 2 where we estimated plant biomass and sampled the soil properties. To select representative farms with traditional and constant grazing use, we interviewed the sheep farmers engaged in grazing in the area. The grazing intensity within the site was determined by the distance to watering places and the number of feces from livestock (Table S1). Grazing herds mainly consisted of Rasa aragonesa sheep.
To test the hypotheses vegetation and soil parameters were selected and determined in those study sites and laboratories. Plant biomass was estimated through nondestructive methods following the proposal by Bonham (2013), in which the average volume of plant species is calculated by multiplying average cover by average height per plot. We calculated the volume of all species in each plot, which were summed up to obtain the total plant biomass in the plot. We ended up with 6 (study areas) Â 3 (grazing level) Â (plots) Â 1 (measurement) replicates (N ¼ 90).
Soil sampling was carried out in two phases. On the one hand, we measured in situ the water infiltration capacity and soil porosity. On the other hand, we collected soil samples (0-10 cm) to measure nitrogen content, microbial biomass, and gypsum content in the soil.
Soil infiltration capacity and porosity were measured within each plot, one under the most abundant plant species and one on bare soil. Infiltration capacity was defined as the value of the one-dimensional cumulative infiltration during 600 s, with the latter calculated from the 1 D cumulative infiltration estimated by applying the measured sorptivity, S, and the hydraulic conductivity, K, to the 1 D Haverkamp et al. (1994) model. S and K were calculated by applying the Latorre et al. (2015) procedure to the corresponding disk infiltrometer infiltration curve measured in the field at 0 cm of soil tension. To this end, we used a 50 mm radius tension infiltrometer disk by Perroux and White (1988). More details about disk infiltrometer measurements can be found in Moret-Fern andez et al. (2021). Total soil porosity was calculated from the total soil water content after saturating an undisturbed soil sample. The soil core (50 mm in diameter and 50 mm in depth) was sampled over the 1-8 cm deep soil layer according to the Grossman and Reinsch (2002) procedure and air-dried for several months to constant weight. Next, the dry soil sample was placed on a sorptivimeter (Moret-Fern andez, Latorre, and Angulo-Mart ınez 2017) and saturated by capillarity up to that the water started to flow by the soil surface. The saturated sample was then weighted, dried at 50 C during 48 h, and weighted again under dry conditions. Total porosity was finally calculated as the difference between the saturated and dry soil weights.
Soil infiltration capacity and porosity are easily understood variables that can be readily related to soil abiotic properties, such as water availability, moisture, or oxygen, which in turn may affect root growth, microbial biomass, or nutrient cycling. A previous study in the same gypseous region by Moret-Fern andez et al. (2021) showed that soil gypsum content had a strong effect on soil hydraulic properties (especially, on soil water infiltration). Thus, we included in the study the gypsum content (57.5% mean between 0.8 and 99.7%) as an independent variable with a potential effect on the plantsoil system under study.
For the determination of soil nitrogen content, microbial biomass, and gypsum content, we extracted three soil samples by digging 15 cm holes in the upper soil profile using a drill, two holes under the two most abundant plant species, and one hole in bare soil. Samples were saved in zip-bags, dried, and stored in a 4 C chamber awaiting lab analyses. In this case, we ended up with 6 (study areas) Â 3 (grazing intensities) Â 5 (plots) Â 3 (analytical) replicates (N ¼ 270) of soil nitrogen content, microbial biomass, and gypsum content in the soil.

Laboratory analysis
The soil nitrogen content was determined by the Dumas method in a Variomax C:N elemental analyzer. Gypsum content was estimated by thermal gravimetric analysis, using the method by Artieda, Herrero, and Drohan (2006).
Microbial biomass was estimated through phospholipid fatty acids analyses (PLFAs), as follows. Soil samples were lyophilized and aliquots of 2 g were used for lipid extraction. Lipids were extracted with one-phase chloroform-methanol-phosphate buffer solvent. Phospholipids were separated from non-polar lipids and converted to fatty acid methyl esters before analysis, following the methodology described by Buyer and Sasser (2012). The resulting fatty acid methyl esters (FAMEs) were separated by gas chromatography using an Agilent 7890 A GC System (Agilent Technologies, Wilmington, DE, USA) equipped with a 25-m Ultra 2 (5%-phenyl)-methylpolysiloxane column (J&W Scientifc, Folsom, CA, USA) and a flame ionization detector. The identification and quantification of FAMEs were carried out using the PLFAD1 method in Sherlock software version 6.3 from MIDI, Inc (Newark, DE, USA). The internal standard 19:0 phosphatidylcholine (Avanti Polar Lipids, Alabaster, AL, USA) was used to quantify the FAMEs. Microbial biomass was estimated by totaling the contents of all individual PLFAs and reported as nanomoles of PLFAs per gram of soil. We ended up with 6 (study areas) Â 3 (grazing level) Â 5 (plots) Â 3 (analytical) replicates (N ¼ 270) of soil nitrogen content, gypsum content, and microbial biomass.

Statistical analysis
Firstly, we obtained the same number of replicates for each measured parameter. We used the mean of the three measurements of soil nitrogen content and microbial biomass per plot, and the mean of the two measurements of infiltration capacity and porosity of soil per plot as the statistic for posterior statistical analyses.

Preliminary analysis
We analyzed the effect of grazing intensity on vegetation and soil properties, i.e., plant biomass, microbial biomass, nitrogen content, infiltration capacity, soil porosity, and gypsum level. We applied one-way ANOVA analyses to test individual models with grazing intensity as the explanatory variable and each vegetation or soil property as a response variable, using the R package car (Fox and Weisberg 2019) in RStudio v.1.1.456 (RStudio Team 2019). Post-hoc Tukey's tests were applied to these models. Parametric assumptions were controlled and variables were transformed if necessary.

Plant-soil interaction model
A general plant-soil interaction model was constructed including the initially hypothesized interactions among grazing intensity, plant biomass, and soil properties (Figure 1a). Plant biomass was modeled as affected by grazing intensity, soil infiltration capacity, and soil porosity; soil nitrogen content as dependent on plant biomass; microbial biomass as dependent on soil nitrogen content, porosity, and infiltration capacity; infiltration capacity as dependent on grazing intensity and gypsum content; and soil porosity as dependent on grazing intensity and gypsum content.
The plant-soil interaction model was analyzed through Structural Equation Modeling (hereafter, SEM). SEM is a statistical technique that allows a set of relationships between one or more independent variables and one or more dependent variables to be tested (Ullman and Bentler 2012). SEM consists of a concatenated analysis of regression models enabling multivariate hypotheses to be tested where variables may act as explanatory and response simultaneously, and providing information on the strength and direction of each interaction. Analyses were performed using the R package lavaan and the "WLSMV" estimator (Rosseel 2012) in RStudio v.1.1.456. The simplified final model (i.e., containing only the significant relationships among variables) was derived by taking out the most significant interactions according to analysis through SEM. This model maintained the same interactions as the initial one, although the direct interaction between nitrogen content and microbial biomass and that between grazing intensity and soil porosity was not included. The amount of variance explained by each interaction and the overall fit was assessed.

Preliminary results
One-way ANOVA showed significant or marginally significant differences among grazing regimes in plant biomass, nitrogen content, microbial biomass, infiltration capacity, and porosity of soil ( Figure 2  According to the post-hoc analysis, plant biomass was significantly different across grazing regimes (Tables 1 and 2). Differences in soil nitrogen content were significant between moderate and high grazing treatments (p < 0.05), in microbial biomass between medium and high grazing, and moderate and high grazing treatments (both p < 0.01), in infiltration capacity between medium and high grazing (p < 0.01), and in porosity between moderate and medium grazing (p < 0.1).

Plant-soil interaction model
Based on the initial GIM output, grazing intensity, soil infiltration capacity, and soil porosity were kept as predictors for plant biomass. The following predictors were kept:  soil nitrogen content, infiltration capacity, and porosity for microbial biomass; grazing intensity for infiltration capacity; and gypsum level for porosity. Final GIM showed a good fit to the null model baseline (Comparative Fit index, CFI: 0.852; Root MeanSquare Error of Approximation, RMSEA: 0.123, p < 0.05, Figure 1b). Grazing intensity had a significant effect on plant biomass. Soil infiltration capacity and porosity had positive and negative effects on plant biomass, respectively. Plant biomass had a significant positive effect on soil nitrogen content but not the reverse. Soil nitrogen content had a significant positive effect on microbial biomass. Microbial biomass was positively and significantly affected by soil nitrogen content and negatively by soil porosity. Grazing intensity significantly affected soil infiltration capacity. Lastly, the gypsum content affected porosity positively.

Discussion
This study investigated the complex set of interactions among grazing intensity, plant biomass, microbial biomass, nitrogen content, porosity, and infiltration capacity in soil. As expected, grazing mainly controlled the plant-soil interaction model through direct changes in plant biomass and soil infiltration capacity. These effects acted indirectly on soil nitrogen content and microbial biomass. Whereas grazing intensity, infiltration capacity, and porosity were the main drivers of plant biomass, microbial biomass was mainly affected by soil nitrogen and porosity.
The hypothesized reduction in plant biomass caused by browsing led to a significant decrease in the soil nitrogen stock, potentially related to reductions in litter input from plant biomass (Harrison and Bardgett 2004). This reduction in nitrogen stock acted indirectly by reducing microbial biomass, suggesting that microbial growth was limited by the amount of nitrogen. Similarly, the negative effect of grazing on soil infiltration capacity had a direct negative effect on plant biomass. Grazers compact the soil by trampling, thus diminishing the pore volume and infiltration capacity of the soil (King and Hutchinson 1976;Moret-Fern andez et al. 2021). The presence of gypsum can enhance this effect by favoring the formation of crust over disturbed areas, as other authors have shown (Shainberg et al. 1989; Moret-Fern andez et al. 2021). These facts limit water infiltration (Eldridge et al. 2020) and thus the availability of water for plants. Whereas soil nitrogen content was strongly dependent on the above-ground plant biomass, plant biomass was strongly dependent on the soil infiltration capacity. These results confirmed our expectations, and they were coincident with previous studies in other arid environments. According to Burke et al. (1998), the most limiting resource for vegetation productivity would be water in arid environments, and nitrogen availability under sub-humid climates. Given that precipitation in this region is <400 mm, water should be the primary control of vegetation productivity, as suggested by our results. This pattern has also been reported in other drylands (Jenny, Smettan, and Facklam-Moniak 1990).
It must be noted that soil infiltration capacity was maximum at medium grazing intensity, which can be related to grazers trampling and crushing the upper soil crust, which would hamper sealing by biological crusts, thereby increasing soil water infiltration (Eldridge et al. 2020;Xiao et al. 2019). This suggests that moderate grazing can partly and indirectly compensate defoliation of plant biomass by enhancing infiltration capacity and thus water availability for plant growth. On the other hand, intensive grazing may have the opposite effect on soil water infiltration, since intensive grazer trampling could favor the formation of compact and impermeable surface crusts that may significantly reduce water infiltration (Wen, He, and Zhang 2016).
Unexpectedly, soil porosity was negatively related to plant biomass and microbial biomass. We expected a positive effect from porosity on plant growth, since a larger pore size could be expected to facilitate root growth, given the smaller resistance to such growth. Similarly, we expected a larger pore size to facilitate fungal expansion in soil (Kravchenko et al. 2011). Other authors found that lower total porosity could prevent water loss by leaching and increase the water retention of gypseous soils (Monson et al. 1992), promoting both root and microbial growth.
Plant and microbial biomass only were indirectly related in the model via nitrogen stock, and thus, we could not verify the hypothesized mutual and direct relationship. However, an indirect facilitation from plants to microbes via nutrient stocks was found. Microorganisms and plant roots can also develop direct and symbiotic relationships to break down the SOM pool when it is recalcitrant (Murphy et al. 2017), which increases nitrogen fixation (Raza et al. 2020) and protects plant roots from pathogens (Shafique et al. 2015).
In conclusion, we found that knowledge of complex direct and indirect plant-soil interactions in drylands and the influence of grazers on the interactions is crucial for the correct management of these areas. In this study, we found that plant-soil interactions were affected by livestock grazing directly through plant biomass and soil infiltration capacity, which indirectly affects the whole plant-soil system. In particular, our results suggest that actions taken to promote soil infiltration and reduce compaction caused by livestock trampling may help in the efficient management of degraded semiarid gypseous rangelands that otherwise would exacerbate the impact of inappropriate management. Moderate grazing regimes are suitable for maintaining the health of rangelands, but grazing pressure needs to be precisely quantified considering both vegetation and soil properties.