Response of arbuscular mycorrhizal fungal community in soil and roots to grazing differs in a wetland on the Qinghai-Tibet plateau

Study site and sampling

The study was carried out in the centre of Zoige Swamp in the Zoige National Nature Reserve on the Qinghai-Tibet plateau (33° 25′−34° 80′N, 102° 29′−102° 59′E, 16,671 ha, 3,365 m above sea level). The site has a plateau cold temperate humid monsoon climate, with a mean annual temperature (MAT) of 1.1 °C, and a mean annual precipitation (MAP) of 660 mm (Wang, Bao & Yan, 2002). The site begins to freeze in late September and is completely thawed in mid-May (Wang, Bao & Yan, 2002). The abundant plant species are Blysmus sinocompressus, Potentilla anserina, Carex enervis, Caltha scaposa, Elymus nutans and Leontopodium wilsonii in the site (Wang, Bao & Yan, 2002).

The Sichuan Zoige Wetland National Nature Reserve Authority approved the collection of soil and root samples in the Zoige Wetland National Nature Reserve. We established 20 plots (each 1 m ×1 m), > 20 m away from each other, in non-grazing (natural grass) and grazing area, respectively (Fig. S1). The average species number and height of vegetation were 7.8 ± 0.495 (mean ± SE) and ca. 31 cm in the non-grazing plots and 5.4 ± 0.255 and ca. seven cm in the grazing plots. This site was mainly grazed by yaks from June to September, as the growing season is from May to August. The grazing intensity (ca. 1.8 yak/ha) in this study site was described as moderate, as previous study showed that there are three different grazing intensities by yaks (light: 1.2 yaks/ha, moderate: 2.0 yaks/ha, and heavy: 2.9 yaks/ha) in alpine meadow on the eastern Tibetan Plateau (Gao et al., 2007). In July 2018, with the best vegetation growth stage, we randomly collected five soil cores (three cm in diameter; 15 cm in depth; ca. 300 g) and mixed into one composite sample from each plot. A total of 40 samples were obtained, packed in an ice box and transported to our laboratory. Soil samples were sieved (1-mm sieve) to remove debris and roots. Subsoil samples were kept at −80 °C until the extraction of fungal hyphae and DNA, and the remaining subsoil samples were air dried and kept at 10 °C until the analysis of AM fungal spore density, GRSP content and soil properties. We manually collected the mixed roots (< two mm in diameter) from each sieved root sample, washed with sterilized deionized water and kept at −80 °C until DNA extraction.

Soil property analysis

We weighed a certain amount of soil from each sample before and after drying for 24 h in an oven at 105 °C, then calculated the percentage of soil moisture. Soil pH was measured at a ratio of 1:2.5 (w/v, soil: water) with a glass electrode (Thermo Orion T20, Columbia, USA). Soil total nitrogen (N) and C were determined by CHNOS Elemental Analyser (Vario EL III Elementar Analysensysteme GmbH, Germany). Soil total phosphorus (P) was extracted using the HClO₄-H₂SO₄ digestion method and determined with a spectrophotometer (UV-2550, Shimadzu, Japan).

EE-GRSP and T-GRSP

Total GRSP (T-GRSP) and easily extracted GRSP (EE-GRSP) were measured according to the method of Janos, Garamszegi & Beltran (2008). We extracted EE-GRSP from 0.1 g air dried soil using sodium citrate buffer (8 mL, 0.02 M, pH 7.0) at 121 °C for 90 min in an autoclave (Yamato SQ810C, China). We repeatedly extracted T-GRSP from 0.1 g air dried soil using sodium citrate buffer (8 mL, 0.05 M, pH 8.0) at 121 °C for 90 min until no obvious color in the supernatant was observed. Supernatants were separated by centrifugation at 6000 g for 15 min to remove the soil particles and saved in a plastic tube (4 °C). Then 0.5 mL of supernatant of EE-GRSP and T-GRSP was stained with 5 mL of Coomassie Brilliant Blue G-250 and was read in a micro-plate reader (Biotek Synergy H4, Winooski, VT, USA) at 595 nm. The bovine serum albumin was used as a standard solution with Coomassie Brilliant Blue method and a standard curve was drawn to determine the content of EE-GRSP and T-GRSP.

AM fungal extra radical hyphal density and spore density

We extracted fungal hyphae from soil according to the membrane filter method (Rillig, Field & Allen, 1999). In total, 4.0 g of frozen soil from each sample was mixed with 12 mL sodium hexametaphosphate (35 g L⁻¹) and 100 mL distilled deionized water in a flask, and then blended for 30 s, settled for 30 min and sieved (38-µm sieve). The fungal hyphae on the sieve were washed into a flask with 200 mL distilled water, and then 2 mL aliquot was filtered through a 25-µm Millipore filter. The fungal hyphae on the filter were stained with 1% acid fuchsine and distinguished into AM and non-AM fungi on the basis of morphological characteristics and staining color (Miller, Jastrow & Reinhardt, 1995). We measured the hyphal length of AM fungi according to the grid-line intersect method (Tennant, 1975). We extracted AM fungal spores from 20 g air dried soil from each sample according to the wet-sieving and decanting method (Daniels & Skipper, 1982) and counted the spore numbers under 40 × magnification (Nikon 80i, Japan).

DNA extraction, PCR and Illumina Miseq sequencing

We extracted DNA from 0.2 g frozen roots and soil using the PowerSoil^® DNA isolation kit (MOBIO Laboratories, Inc., Carlsbad, USA) in accordance with the manufacturer’s instructions, and measured the DNA concentration using a NanoDrop 1000 Spectrophotometer (Thermo Scientific, Wilmington, USA). We amplified the fungal 18S rDNA region using a two-step PCR procedure. The first PCR using primers AML2 (Lee, Lee & Young, 2008) and GeoA2 (Schwarzott & Schüßler, 2001) was conducted in a final 25 µL reaction mixture, including ca. 10 ng of template DNA, 0.75 µM of each primer, 250 µM of each dNTP, 0.5 U KOD-plus-Neo polymerase (Toyobo, Tokyo, Japan), 1.5 mM MgSO₄, and 2.5 µL 10 × buffer. The thermal cycling conditions were performed as follows: an initial denaturation at 95 °C for 5 min, 30 cycles for denaturation at 94 °C for 1 min, annealing at 58 °C for 50 s and extension at 68 °C for 1 min, and a final extension at 68 °C for 10 min. The products of the first amplification were diluted 100 times, and 1 µL of the diluted DNA template was used for the second amplification. The thermal cycling conditions for the second amplification were the same as first amplification, except that the primers NS31 (Simon, Lalonde & Bruns, 1992) and AMDGR (Sato et al., 2005) linked with 12-base barcode sequences were used. The size of amplified fragment was about 300 base pairs (bp). We purified the PCR products using a PCR Product Gel Purification Kit (Omega Bio-Tek, USA), and pooled the purified PCR products with the same amount (100 ng) from each sample and adjusted the concentration to 10 ng µL⁻¹. We constructed a sequencing library by addition of an Illumina sequencing adaptor (5′-GATCGGAAGAGCACACGTCTGAACTCCAGTCACATCACGATCT- CGTATGCCGTCTTCTGCTTG-3′) to the products using the Illumina TruSeq DNA PCR-Free LT Library Prep Kit (Illumina, CA, USA) according to the manufacturer’s instructions. We sequenced the library by an Illumina MiSeq PE 250 platform using the paired-end (2 × 250 bp) option in the Chengdu Institute of Biology, Chinese Academy of Sciences, China.

Bioinformatics analysis

We filtered the raw sequences using Quantitative Insights into Microbial Ecology (QIIME) v.1.7.0 (Caporaso et al., 2010) to eliminate low-quality sequences, such as read length < 200 bp, no valid primer sequence or barcode sequence, containing ambiguous bases, or an average quality score < 20. We checked and deleted the potential chimeras against the MaarjAM database (Öpik et al., 2010) using the ‘chimera.uchime’ command in Mothur version 1.31.2 (Schloss et al., 2009). High quality sequences were subjected to de-replication and de-singleton, and then clustered into operational taxonomic units (OTUs) at a 97% sequence similarity level using the cluster_otus command in USEARCH v8.0 (Edgar, 2013). Using a basic local alignment search tool (BLAST) (Altschul et al., 1990), we selected the most abundant sequence of each OTU and searched against the MaarjAM database and National Center for Biotechnology Information (NCBI) nt database. We identified OTUs as the AM fungi based on the closest BLAST hit annotated as ‘Glomeromycotina’ and E values < e⁻⁵⁰. Furthermore, we normalized the sequence number of each sample to the smallest sample size using the ‘sub.sample’ command in Mothur. We have submitted the representative sequence of each AM fungal OTU to the European Molecular Biology Laboratory (EMBL) database (accession no. LR736402-LR736557). The identified AM fungi are shown in Table S1 .

Statistical analysis

We conducted all statistical analyses in R version 3.3.2 (R Development Core Team, 2017). Tukey’s honestly significant difference (HSD) test or Conover’s test was used to examine the significant difference of soil moisture, pH, total N, total C and total P in the grazing and non-grazing plots at P < 0.05. Generalized linear model (GLM) with Poisson error structure and log link function was conducted to evaluate the effect of grazing on AM fungal spore density, extra radical hyphal density and T-GRSP, as these data did not meet normal distribution, while GLM with Gaussian error structure and identity link function was conducted to evaluate the effect of grazing on EE-GRSP, and then Conover’s test was used to examine the significant difference between grazing and non-grazing treatments at P < 0.05 using the post-hoc.kruskal.conover.test function in the PMCMR package (Pohlert, 2014). Meanwhile, GLM with Gaussian error structure and identity link function was used to evaluate the effect of grazing and sample type nested in grazing on the AM fungal OTU richness, and GLM with Gamma error structure and inverse link function was used to evaluate the effect of grazing and sample type nested in grazing on the relative abundance of abundant OTUs (relative abundance > 1%) and orders of AM fungi, as these data did not meet normal distribution, and then Conover’s test was conducted for comparisons between grazing and non-grazing treatments in soil and roots at P <  0.05.

The distance matrices of AM fungal community composition (Hellinger-transformed OTU read data) in roots and soil were established by the Bray–Curtis method (Clarke, Somerfield & Chapman, 2006). Nested permutational multivariate analysis of variance (PerMANOVA) was conducted to examine the effect of grazing and sample type nested within grazing on AM fungal community composition, using the ‘adonis’ function in the vegan with 999 permutations (Oksanen et al., 2013). Besides, PerMANOVA was conducted to examine the effect of grazing on AM fungal community composition in soil and roots, respectively. Redundancy analysis (RDA) was conducted to reveal the significant correlation of AM fungal community composition and soil variables using the Monte Carlo permutation test with 999 permutations.

Soil properties

Soil pH was significantly lower in grazing treatment than in non-grazing treatment (Table 1). Soil moisture, total N, total C and total P were significantly higher in grazing treatment than in non-grazing treatment (Table 1).

EE-GRSP and T-GRSP contents

The EE-GRSP content was 28.96  ± 3.73 µg g⁻¹ (mean  ± SE) and 25.71  ± 2.26 µg g⁻¹ in grazing and non-grazing treatments, respectively. The T-GRSP content was 96.7  ± 18.82 µg g⁻¹ and 71.87   ± 12.87 µg g⁻¹ in grazing and non-grazing treatments, respectively. GLM showed that grazing significantly influenced EE-GRSP (P = 0.002; Table 2) and T-GRSP (P < 0.001; Table 2). For example, EE-GRSP and T-GRSP contents were significantly lower in non-grazing than in grazing treatments (Figs. 1A and 1B).

AM fungal spore density and extra radical hyphal density

The spore density of AM fungi was 25.89  ± 12.17 g⁻¹ (mean  ± SE) and 15.03 ± 5.88 g⁻¹ in grazing and non-grazing treatments, respectively. The extra radical hyphal density of AM fungi was 4.00  ± 2.51 m g⁻¹ and 3.10  ± 1.56 m g⁻¹ in grazing and non-grazing treatments, respectively. GLM revealed that grazing significantly affected AM fungal spore density (P = 0.001; Table 2) but not extra radical hyphal density (P = 0.130; Table 2). For example, the spore density of AM fungi was significantly lower in non-grazing than in grazing treatments (Fig. 1C). However, AM fungal extra radical hyphal density was not significantly different in non-grazing and grazing treatments (Fig. 1D).

Identification of AM fungi

In total, 3,205,557 high-quality sequences were filtered from 3,335,816 raw sequences and clustered into 882 OTUs at a 97% sequence similarity level. Among 882 OTUs, 156 (2,919,706 sequences) belonged to AM fungi. As the sequence number of AM fungi varied from 20,408 to 48,572 in the 80 samples, the number of sequence was normalized to 20,408. The normalized dataset contained 156 AM fungal OTUs (1,632,640 sequences). Of the 156 AM fungal OTUs obtained, 154 were from soil, 152 from roots, and 150 shared both soil and roots. Among 156 AM fungal OTUs, 153 were detected from more than three samples (frequency ≥ 3.75%) (Fig. S2A). Furthermore, the 21 abundant AM fungal OTUs (relative abundance > 1%) occupied 83.85% of the total sequences (Fig. S2B). Among 156 AM fungal OTUs, 109 were identified to Glomerales (79.52% of sequences), 22 to Diversisporales (10.84%), 21 to Archaeosporales (8.75%), and 4 to Paraglomerales (0.89%). In addition, the rarefaction curves indicated that the sample numbers were sufficient to detect the most AM fungi in this study (Fig. S3).

AM fungal OTU richness

AM fungal OTU richness in grazing and non-grazing treatments was 123.70  ± 2.96 (mean  ± SE) and 122.85  ± 3.01 in soil, and 117.55  ± 2.26 and 118.00  ± 4.38 in roots, respectively. GLM revealed that AM fungal OTU richness was influenced by sample type (root and soil; P < 0.001; Table S2), but not by grazing (P = 0.662; Table S2). For example, AM fungal OTU richness was significantly lower in roots than in soil in both grazing and non-grazing treatments (Fig. 2). However, AM fungal OTU richness was not significantly different between grazing and non-grazing treatments in roots and soil, respectively (Fig. 2).

AM fungal community

GLM revealed that grazing had significant effect on the relative abundance of abundant AM fungal OTU12 and OTU25 (Glomerales), and sample type had significant effect on the relative abundance of abundant AM fungal OTU4, OTU5, OTU7, OTU12, OTU14, OTU18, OTU25 and OTU141 (Glomerales), OTU8 and OTU17 (Diversisporales) and OTU23 (Archaeosporales) (Fig. 3; Tables S3 ; S4).

GLM revealed that sample type significantly influenced the relative abundance of Glomerales, Diversisporales and Archaeosporales, and grazing significantly affected the relative abundance of Diversisporales (Fig. 4; Table S5). The relative abundance of Glomerales was significantly lower in soil than in roots; by contrast, the relative abundance of Diversisporales and Archaeosporales was significantly lower in roots than in soil, regardless of non-grazing and grazing treatments (Fig. 4; Table S5). Besides, the relative abundance of Diversisporales was significantly lower in grazing treatment than in non-grazing treatment (Fig. 4; Table S5).

The PerMANOVA demonstrated that the community composition of AM fungi was significantly influenced by sample type (soil and root; F = 7.2836, R² = 0.157, P = 0.001; Table S6) and grazing (F = 2.339, R² = 0.025, P = 0.012; Table S6). Furthermore, the community composition of AM fungi was significantly influenced by grazing in soil (F = 2.639, R² = 0.055, P = 0.001; Table S7), but not in roots (F = 0.998, R² = 0.025, P = 0.419; Table S8). Furthermore, RDA showed that the community composition of AM fungi in soil and roots was significantly correlated with soil pH, moisture, total C, total N and total P (Fig. 5).