Responses of fungal communities at different soil depths to grazing intensity in a desert steppe

Experimental area

The experimental area is located in Zhurihe Town, Sunite Right Banner, Xilin Gol League, Inner Mongolia, China (112°47′16.9″E, 42°16′26.2″N). The area belongs to the desert steppe ecoregion, with a soil type of light chestnut calcisols. The thickness of the humus layer ranges from 5 to 10 cm. The vegetation in the experimental area is dominated by the species Stipa breviflora as foundation species, with Allium polyrhizum and Cleistogenes songorica as dominant species. Other common species include Convolvulus ammannii, Kochia prostrata, Allium tenuissimum, Heteropappus altaicus, Cleistogenes squarrosa, Caragana stenophylla, and Carex duriuscula. The annual precipitation in the experimental area ranges from 170 mm to 190 mm, with an average annual temperature of 5.8 °C, classifying it as a mesothermal climate.

Experimental design

A total of five treatments were set up in this experiment, with three replicates for each treatment, resulting in 15 experimental grazing plots. Samples were taken at two different locations and two different depths at each test site. A completely randomized block design was employed. The total area of the grazing experimental site was 39.26 hectares. The five grazing treatments were: no grazing (NG), light grazing (LG), moderate grazing (MG), heavy grazing (HG), and over grazing (OG). Based on the local government’s guidance documents, as well as the actual grazing practices observed among the herdsmen, the grazing intensities for the various treatment areas were determined. Specifically, the NG treatment had 0 sheep, while the LG, MG, HG, and OG treatments had four, five, six and seven sheep, respectively. These intensities translate into stocking rates of 0 sheep/ha, 1.54 sheep/ha, 1.92 sheep/ha, 2.31 sheep/ha, and 2.69 sheep/ha, respectively.

Sample collection and index determination

The experimental site was established in 2022, and the grazing period lasts for six months, starting in mid-May and ending in mid-October each year. In May 2023, three different sampling points were randomly selected in each experimental plot. We used a soil auger (with a diameter of 10 cm) to collect soil samples from both the topsoil (0–20 cm) and the subsoil (20–40 cm). Based on previous studies, we discovered that the soil texture undergoes significant changes at a depth of 20 centimeters. Therefore, we selected this depth as the boundary for our study. The three soil samples from the same plot were mixed together to form a single sample, and then passed through a 2-mm sieve to remove plant roots and gravel. The mixed samples were divided into two portions. One portion was placed in sterile centrifuge tubes inside an ice box and transported back to the laboratory for storage at −20 °C for the determination of soil microbial community structural characteristics. The other portion was placed in a self-sealing bag for the measurement of soil physicochemical properties.

Index determination of soil physical and chemical properties

To assess soil physicochemical properties, water and soil samples were mixed at a ratio of 2.5:1, and the supernatant of the mixed solution was used to measure the pH using a pH meter. The soil organic matter (SOM) was calculated based on the content of soil organic carbon (SOC), which was determined by oxidation with potassium dichromate and external heating (Nelson, Sommers & Sparks, 1996). To assess dissolved organic carbon (DOC) content, water and soil samples were mixed at a ratio of 5:1, filtered through a 0.045 µm membrane, and measured using a TOC-L instrument. After extracting the soil sample with potassium chloride, soil nitrate nitrogen (NO${}_3^{-}$-N) and ammonium nitrogen (NH⁴⁺-N) were determined using a flow analyzer. Soil total nitrogen (TN) was determined using the semi-micro Kjeldahl method. β-1,4-glucosidase (βG), N-acetyl-β-D-glucosidase (NAG), and alkaline phosphatase (AKP) were determined using the microplate fluorescence method (Marx, Wood & Jarvis, 2001). SWC was determined by gravimetric method. First measure the weight m0 of the drying aluminum box, put a certain amount of fresh soil into the aluminum box, weigh the weight m1 of the fresh soil and the aluminum box, and measure the weight m2 of the drying soil and the aluminum box after drying.

DNA extraction and sequencing of soil fungi

The genomic DNA of the sample was extracted using the CTAB method (Pahlich & Gerlitz, 1980), and then the purity and concentration of the extracted DNA were detected using agarose gel electrophoresis. Subsequently, an appropriate amount of DNA was taken into a centrifuge tube and diluted to 1ng/µl with sterile water. Using the diluted genomic DNA as a template, PCR was performed with barcoded specific primers targeting the amplification region and Phusion^® High-Fidelity PCR Master Mix (New England Biolabs, Ipswich, MA, USA) to ensure amplification efficiency and accuracy. PCR amplification of the samples was performed using ITS1 region primers (F: CTTGGTCATTTAGAGGAAGTAA; R: GCTGCGTTCTTCATCGATGC). PCR products were detected using 2% agarose gel electrophoresis. The PCR products were mixed in equal proportions, and then Qiagen Gel Extraction Kit (Qiagen, Hilden, Germany) was used to purify the mixed PCR products. Library preparation was carried out using the NEBNext^® Ultra™ II DNA Library Prep Kit. The library quality was evaluated on the Qubit@ 2.0 Fluorometer (Thermo Fisher Scientific, Waltham, MA, USA) and Agilent Bioanalyzer 2100 system. Once the libraries met the required standards, they were sequenced on the NovaSeq 6000 platform. Paired-end reads were merged using FLASH (Version 1.2.11), and the splicing sequences were called Raw Tags. Quality filtering on the raw tags were performed using the fastp (Version 0.20.0) software to obtain high-quality Clean Tags. The Clean Tags were compared with the reference database (Unite database https://unite.ut.ee/ for ITS) using Vsearch (Version 2.15.0) to detect the chimera sequences, and then the chimera sequences were removed to obtain the Effective Tags (Haas et al., 2011). The effective sequencing data obtained was denoised using DADA2, and sequences with an abundance less than 5 were filtered out (Li et al., 2020), resulting in the final amplicon sequence variants (ASVs). The average number of reads per sample exceeded 35,000. Species annotation was performed using QIIME2 software.

Data analysis

We used Microsoft Excel 2021 and SPSS to initially organize the data, and used R 4.2.2 to analyze and plot the data. To determine the effects of different grazing intensities and soil depth on soil physicochemical indicators (SOM, TN, ammonium nitrogen, nitrate nitrogen, DOC, pH, and SWC), soil enzyme activities (βG, AKP, and NAG), and soil fungal communities (abundance and α-diversity), a one-way analysis of variance (ANOVA, Tukey test) and T test was conducted. By analyzing these specific soil parameters under different grazing intensities, researchers can gain a better understanding of how these management practices impact soil health and fungal ecology.

Microbial diversity metrics, including total Observed_otus, Chao 1, Shannon diversity, Simpson diversity, Pielou_e, and Goods_coverage, were calculated using the vegan package in R software and visualized using the ggplot2 package. Chao 1 and Observed_otus indices are used to describe species richness, while Shannon and Simpson indices are measures of species diversity. The Pielou_e index is used to describe species evenness. Multivariate PERMANOVAs analysis was conducted using the vegan package in R software. Principal coordinates analysis (PCoA) based on Bray-Curtis distances was performed, and the results were visualized using the ggplot2 package. PCoA could show the differences in fungal communities, and PERMANOVAs analysis proved the statistical characteristics of the differences. It demonstrates the relationship between microbial communities at different depths and under varying grazing intensities. Spearman correlation between soil physicochemical indicators, soil enzyme activity, and soil fungi was calculated using the psych package in R software and visualized accordingly. Heatmaps illustrate the relationship between soil properties and fungal phyla at two soil depths. Multivariate regression trees (MRT) analysis, utilizing the mvpart package in R software, was employed to assess the most significant abiotic factors influencing fungal α-diversity (total observed_otus, Chao 1, Shannon diversity, Simpson diversity, Pielou_e) and community composition (fungal phyla) (De’Ath, 2002).

Soil characteristics

Grazing affects soil characteristics, but the specific nutrients and the magnitude of differences depend on soil depth. As shown in Table 1, in topsoil, LG had 10.8% higher SOM content than NG, MG had 20.3% lower DOC than NG, LG and HG had higher SWC than NG (35.4% and 50.8% higher, respectively), while LG, MG, and OG had lower AKP content than NG (29.9%, 56.5%, and 30.8% lower, respectively). OG had a 2.8% lower TN content than NG. MG, HG, and OG had lower βG contents than NG (61.1%, 46.3%, and 61.3% lower, respectively), while all treatments including grazing had lower NAG contents than NG (17.9%, 33.6%, 40.5%, and 74.6% lower, respectively. In the subsoil, HG had higher SOM and TN than NG (12.3% and 10.0% higher, respectively). LG, HG and OG had lower AKP contents in the soil than NG (15.8%, 42.9%, and 45.2% lower, respectively), and the four treatments involving grazing had lower βG contents than NG (30.0%, 39.7%, 24.9%, and 39.1% lower, respectively, while the NAG contents were 46.1%, 41.3%, 50.7%, and 60.3% lower than NG, respectively. OG had higher nitrate nitrogen that NG (308.2% higher in the topsoil and 466.2% higher in the subsoil). The SWC of the subsoil was higher than that of the topsoil in NG, LG, MG, and HG treatments, but NAG showed the opposite trend. Except for MG, AKP had the same trend as NAG. The content of DOC in surface soil of HG was significantly higher than that of subsurface soil. The βG content in surface soil of control treatment and grazing treatment was higher than that of subsurface soil. Soil depth reduces the effect of grazing on soil physicochemistry.

Fungal α-diversity

Like for the soil characteristics, the effects of grazing on fungal α diversity were different at different depths. In the topsoil, the effect of grazing on fungal α-diversity was not significant (Fig. 1). In the subsoil, the fungal α diversity indices in all treatments involving grazing were lower than those in the plots without grazing, and the Observed_otus index, and Chao1 index in LG, MG and HG were significantly lower than those in the NG treatment (Fig. 1). There were no differences among treatments involving grazing. The Pielou index and Simpson index were also lower, but not significantly (Fig. 1). The Shannon indices of LG and OG were significantly lower than the index for NG (Fig. 1). For both topsoil and subsoil, all grazing treatments reduced α diversity indices of the subsoil fungi. Grazing also changed the effect of depth on fungal α diversity. In the NG treatment, the Observed_otus index and Chao1 index in the topsoil were lower than those in the subsoil (Fig. 1). However, there were no significant differences in fungal α diversity at different depths under the different treatments with grazing.

Fungal β-diversity

Grazing had a large effect on the relative abundance of fungi, and the effects of grazing on fungal relative abundance were different at different depths. The PCoA sorting results showed that there were significant differences in fungal community composition (beta diversity) in the treatments involving grazing than in the NG. There were also significant differences in fungal community composition between light, moderate, and heavy grazing (Fig. 2). There were no differences in fungal community composition between the two soil depths, similarly, bifactor analysis in grazing and soil depth did not find significant differences (Fig. 2).

The dominant phyla were Ascomycota, Basidiomycota, Mortierellomycota and Glomeromycota. In the topsoil, the relative abundance of dominant fungi was not significantly different. In the subsoil, the relative abundance of Mortierellomycota was lower in LG, HG, and OG than in the NG by 82.7%, 61.9%, and 47.7%, lower, respectively; (Fig. 2B). Glomeromycota was 51.5% higher in HG than NG (Fig. 2B). Compared with topsoil, Glomeromycota in subsoil was reduced by 20.4% (Table S3). However, the relative abundance of Chytridiomycota and Calcarisporiellomycota in subsoil was 218% and 337%, respectively, of that in topsoil (Table S3).

At the genus level, in the topsoil, compared with NG, HG and OG make the relative abundance of Coprinopsis decreased by 99.5% and 98.8%, respectively (Table S1). LG,MG,HG and OG significantly decreased the relative abundance of Clitocella (Table S1). LG and MG significantly decreased the relative abundance of Limonomyces (Table S1). HG significantly increased the relative abundance of Ceratobasidium (Table S1). OG caused a significant decline in the relative abundance of Cercophora (Table S1).

In the subsoil, compared with NG, LG and HG significantly decreased the relative abundance of Mortierella and Cercophora (Table S1). HG significantly increased the relative abundance of Glomus and Ceratobasidium (Table S1). LG,MG,HG,OG significantly decreased the relative abundance of Preussia and Clitocella (Table S1). OG significantly decreased the relative abundance of Coprinopsis (Table S1). LG significantly decreased the relative abundance of Darksidea and Limonomyces (Table S1). However, there was no significant difference in the relative abundance of fungi genera at different soil depths.

The linkage between soil characteristics and fungal communities

The correlation between fungi and soil environmental factors varied with soil depth, as heatmap analysis revealed significant relationships between soil physicochemical properties and fungi, with distinct patterns observed in the topsoil and subsoil. The topsoil and subsoil showed different patterns. In the topsoil, Chytridiomycota were significantly correlated with DOC, βG, AKP, and SWC, Basidiomycota were significantly positively correlated with AKP, while Kickxellomycota were significantly positively correlated with pH (Fig. 3). Calcarisporiellomycota were significantly positively correlated with SWC (Fig. 3). In the subsoil, pH significantly correlated with Basidiomycota, Glomeromycota, and Calcarisporiellomycota, while AKP had a significant negative correlation on Glomeromycota (Fig. 3). Rozellomycota had a significant negative correlation with TN and SOM (Fig. 3). Basidiomycota had a significant positive correlation with NAG (Fig. 3). Calcarisporiellomycota were significantly positively correlated with SWC (Fig. 3).

Through multiple regression tree analysis, the main driving factor regulating fungal community changes appeared to be SWC in the topsoil and ammonium nitrogen and nitrate nitrogen in the subsoil (Figs. 4A and 4B). For changes in fungal α-diversity, TN was the main factor influencing the topsoil, while nitrate nitrogen had a significant impact on the subsoil (Figs. 4C and 4D).

Impact of grazing on soil physiochemical properties

Our results demonstrate that the impact of different grazing intensities on topsoil and subsoil is not entirely similar. In the topsoil, soil organic matter (SOM) content was affected by grazing, with first an increasing and then a decreasing trend with an increase in intensity, reaching its maximum under light grazing. This finding is consistent with previous research (Song et al., 2023), which indicated that moderate grazing can enhance SOM content. Light grazing can have a positive effect on soil organic matter content through production of manure, increase decomposition in the plant litter, and activation of nutrients by root exudates (Inagaki et al., 2023; Mikutta et al., 2019; Sokol, Sanderman & Bradford, 2019). Our results show grazing reduce the plant litter (Table S2). Grazing livestock can impact plants through various behaviors such as feeding and rubbing, which causes changes in plant litter. Additionally, trampling by livestock accelerates the decomposition of litter in the soil. However, as grazing intensity increases, its effect on litter decomposition gradually decreases (Su et al., 2022; Wei et al., 2023).

Previous studies revealed that the impact of grazing intensity on soil total nitrogen content is nonlinear (Dai et al., 2022). In grasslands, the total nitrogen content initially increases with grazing intensity and then decreases, reaching a maximum under light grazing (Shi et al., 2023; Cao et al., 2024). This is not completely consistent with our results, which may be related to the type of grassland or the climate environment, but the overall trend is consistent. Livestock manure and plant litter are two main sources of soil nitrogen (Wang et al., 2020), the decomposition rate of manure and litter in the soil decreases as grazing intensity increases. Nitrate nitrogen in the soil mainly originates from the urine of livestock. Previous studies showed that as grazing intensity increases, the content of soil nitrate continuously increases (Li et al., 2023a; Li et al., 2023b; Zhang et al., 2023a; Zhang et al., 2023b; Lv et al., 2024). This is consistent with our results. An increase in grazing intensity leads to a decrease in plant biomass (Table S2), leading to decreased absorption of nitrate by plants (Kant, 2018), Additionally, grazing can also reduce the loss of nitrate from the soil (Eriksen et al., 2015), resulting in accumulation of nitrate in the soil.

Compared with NG, soil DOC content became lower under LG, MG, and OG conditions, but higher under HG. Soil DOC originates from plant litter, and the effects of grazing on DOC are similar to the above effects of grazing on SOM. The influence of soil moisture is more pronounced during periods of low precipitation (Luo et al., 2009). Our results show Grazing increased soil water content, reaching a maximum under MG. Grazing increases surface litter, but livestock trampling accelerates the decomposition rate of the litter. Previous studies found that the increase in litter can reduce soil water evaporation, thereby increasing SWC (Zhao et al., 2023). Due to the effect of litter and water evaporation, soil water content reaches its maximum under MG.

In the subsoil, the trends of nitrate, DOC, SWC, SOM, and TN were consistent with those in the topsoil. However, the content of SOM and TN reached its maximum under HG, rather than under LG. This is because the SOM in surface soil layers mainly originates from root exudates and not from plant litter (Jobbágy & Jackson, 2000), and the impact of livestock on soil decomposition capacity, jointly determine this result. Plant growth requires nitrogen (N) absorption through their roots, and moderate grazing promotes plant growth (Song et al., 2022). Therefore, under LG and MG, TN content decreases due to the utilization of N by plants. But when the number of plants falls to a certain number, the reproduction ability declines. Therefore, under HG and OG, excessive herbivory by livestock leads to a decrease in the number of plants, reducing the absorption of N and resulting in an increase in TN content (Table S2). Compared to the topsoil, the subsoil generally has lower nutrient content, with the exception of SWC. This is due to our experimental site being located in a desert steppe, where the water-holding capacity of subsoil is stronger than that of the topsoil. For nitrate, DOC, and SWC, the influence of grazing is mitigated by soil depth, as livestock activities have a greater impact on topsoil than on subsoil.

Impact of grazing and soil depth on fungi

Fungi are heterotrophic organisms that primarily feed on decaying substrates. Grazing has an impact on both the quality and quantity of litter, which can easily alter fungal diversity in grasslands (Rong et al., 2022). Grazing treatments can also affect fungi by altering the interactions between soil fungi and their habitats (Guitian & Bardgett, 2000; Schipper et al., 2014; Zhang et al., 2020). Our results show that various α-diversity indices of soil fungi under treatments with grazing are lower than those under the NG treatment, which is inconsistent with previous findings (Wang & Tang, 2019). On the one hand, grazing alters the soil environment, leading to changes in fungal diversity. For example, grazing altered SOC thus affecting soil microbial communities (Zhang et al., 2024). On the other hand, grazing treatment modifies plant species composition, resulting in changes in the fungal community recruited through root exudates (Xu et al., 2022). Changes in plant species composition, such as tomato and onion intercropping, increased the colonization of specific Bacillus spp and established specific rhizosphere microbial communities compared with tomato cultivation alone (Zhou et al., 2023). The various α-diversity indices in subsoil were higher than those in topsoil. These results are consistent with previous reports (Wang et al., 2023). This is because, compared to topsoil, fungi in subsoil are stronger influenced by root exudates. During plant growth, root exudates exert a recruitment effect (Zhou et al., 2023), resulting in a higher α-diversity of fungi in subsoil than in topsoil.

In our study, the primary dominant fungal phyla in all soils were Ascomycota, followed by Basidiomycota, Glomeromycota, and others. Our experimental site belongs to a desert steppe, where fungi require stronger decomposition abilities to obtain organic matter than in other environments. Ascomycota are the primary decomposer among soil fungi (Inagaki et al., 2023). Ascomycota is the most dominant fungal phylum, with a relative abundance exceeding 40% in all treatments. Furthermore, Ascomycota are more abundant in subsoil than in topsoil. The reason for this is that topsoil is rich in humus, which is easy to decompose. Compared to topsoil, subsoil contains a greater variety of recalcitrant substances. Therefore, more Ascomycota are required in subsurface soil to decompose these substances (Manici et al., 2024). Since subsoil has less soil OM than topsoil, oligotrophic fungi such as Basidiomycota are more abundant in topsoil, while copiotrophic fungi such as Ascomycota are more abundant in subsoil (Yang et al., 2023). Grazing alters the relative abundance of fungi. For instance, under LG, HG, and OG treatments, the relative abundance of Mortierellomycota was significantly lower due to changes in the soil environment caused by grazing (Ahonen et al., 2024). Additionally, with increasing grazing intensity, certain fungal genera belonging to the same phylum exhibited opposite changes. For example, within Basidiomycota, Ceratobasidium and Coprinopsis displayed opposite trends. This may be due to the fact that these two fungal genera within the Basidiomycota phylum play similar ecological roles in the soil (Snajdr et al., 2011), but they prefer different soil environments for growth.

Regulating effects of soil properties on the response of fungi to grazing intensity

Multivariate regression tree (MRT) analysis revealed that soil total nitrogen (TN) was the primary factor driving changes in fungal α-diversity in the topsoil, while nitrate nitrogen was the main factor in the subsoil. Previous studies also demonstrated that nitrogen elements in the soil significantly influence changes in soil fungi (Chen et al., 2021). In the subsoil, βG showed a positive correlation with the Chao1 and Observed_otus indices (Table S1). βG is a key enzyme for microbial carbon acquisition (Sinsabaugh et al., 2008). Compared with topsoil, fungi in the subsoil require more carbon to maintain their biological functions (Huang et al., 2023). The MRT analysis revealed that the main driving factor regulating fungal community changes was SWC in the topsoil, while it was ammonium nitrogen and nitrate nitrogen in the subsoil. For the topsoil, previous studies have shown that soil microorganisms are very sensitive to nutrients and soil water, and even small changes in water content can lead to changes in the dominant members of the soil microbial community (Ramírez et al., 2020; Cao et al., 2022). When soil water content decreases, soil aeration increases (Levy-Booth et al., 2016); this may affect the relative abundance of anaerobic and aerobic bacteria in the soil. In the subsoil, nutrients such as nitrogen have a strong influence on the fungal community due to the lack of nutrients.

In our correlation analysis, Glomeromycota in subsoil showed a strong negative correlation with AKP (alkaline phosphatase) content and a positive correlation with pH. Glomeromycota, also known as arbuscular mycorrhizas fungi, can promote plant absorption of phosphorus (Bushra et al., 2024). The decrease in phosphorus (P) content affects the AKP levels in the soil (Gao et al., 2024). Additionally, plants absorb P from the soil in the form of phosphates and secrete various substances to enhance soil nutrient content (Oldroyd & Leyser, 2020), which in turn promotes fungal growth. With increasing soil depth, the correlation of soil characteristics such as DOC, SWC, AKP and βG with Chytridiomycota disappeared, and negative correlations of TN and SOM with Rozellomycota emerged. This may be because changes in the nutrient environment result in altered fungal life strategies (Yang et al., 2023). The change in fungal life strategy liberates the fungus from being previously limited by nutrients or creates new constraints. Nutrient-rich fungi such as Chytridiomycota typically exhibit a positive correlation with soil nutrients when soil nutrient content is low. However, in conditions of abundant soil nutrients, they may be subject to self-limitation. The oligotrophic fungi such as Rozellomycota, on the other hand, show a negative correlation with soil nutrients when soil nutrient content is high.