Alterations to arbuscular mycorrhizal fungal community composition is driven by warming at specific elevations

Site description

The Qinghai-Tibet Plateau is a vast plateau in Central Asia covering most of the Tibet Autonomous Region and Qinghai Province in China. Referred to as the “the roof of the world” occupying 2.5 million square kilometers, it is the highest and biggest plateau of the world. The annual average temperature is −4 °C, with annual precipitation ranges from 100 to 300 mm. In our study, the main vegetation is Kobresia pygmea and the type of soil is alpine meadow soil. The slope for each sampling site is less than 2°. In view of the uniqueness of climatic and topographical characteristics on QTP, this study selected samples between 29°50′36.49″–29°54′26.70″ north latitude and 102°0′42.50″–102°02′9.50″ east longitude on the eastern part of QTP (Table 1).

Experiment design and sample collection

Quadrats of 20 m × 20 m were positioned at four elevations of 3,000, 3,500, 3,800 and 4,170 m on QTP. Each quadrat was divided into 25 of 4 m × 4 m little quadrats. We took three biological repetitions with non-adjacent randomly as control treatment (CK) and OTC warming treatment (OTC) by the way of artificial and simulated warming through open-top chamber, respectively (Gao & Li, 2019; Li et al., 2020a). Compared with other warming methods, it can ensure that the soil is basically undamaged and easily to repeated (Klein, Harte & Zhao, 2004). The top and bottom are hexagonal and open with the side composed of six trapezoid-shaped plexiglass. We carried out a 1-year warming test and all samples were taken in August and September of the next year without rain or snow. Small meteorological observation stations were set up at each altitude to monitor soil temperature and moisture. Instantaneous measurement of soil temperature and soil moisture was performed by fixed-point measurement using HOBO PRO temperature and soil moisture recorder. We selected soils samples randomly with a soil corer with diameter of two cm and depth of 0–20 cm. We mixed three soil cores as a sample and repeated three times in CK and OTC, respectively. Then, separating the root system from the soil and sealing with sterile plastic valve bags, with DNA samples being stored at –20 °C. Field experiments were approved by the Key Laboratory of Mountain Surface Processes and Ecological Regulation, Chinese Academy of Sciences (20160416).

DNA extraction and PCR amplification

Genomic DNA was extracted from soil samples, using the Fast DNA SPIN Kit for Soil (MP Biomedicals LLC, Santa Ana, CA, USA) according to manufacturer’s protocols. The final DNA purification and concentration were determined by Nano Drop 2000 UV-vis spectrophotometer (Thermo Scientific, Wilmington, NC, USA), and DNA quality was checked by 1% agarose gel electrophoresis. The extracted DNA was subjected to nested PCR by thermocycler PCR system (GeneAmp 9700; ABI, Sunnyvale, CA, USA). The first PCR amplification was performed with primers AML1F (5′-ATCAACTTTCGATGGTAGGATAGA-3′) and AML2R (5′-GAACCCAAACACTTTGGTTTCC-3′) by an ABI GeneAmp® 9700 PCR thermocycler (ABI, Sunnyvale, CA, USA) (Lee, Lee & Young, 2008; Li et al., 2019). The PCR reactions were conducted using the following program: 3 min of denaturation at 95 °C, 32 cycles of 30 s at 95 °C, 30 s for annealing at 55 °C, and 45 s for elongation at 72 °C, and a final extension at 72 °C for 10 min. PCR reactions were performed in triplicate 20 μL mixture containing 4 μL of 5 × FastPfu Buffer, 2 μL of 2.5 mM dNTPs, 0.8 μL of each primer (5 μM), 0.4 μL of FastPfu Polymerase and 10 ng of template DNA. The second PCR amplication used identical reaction conditions described above with the primers AMDGR (5’-CCCAACTATCCCTATTAATCAT-3’) and AMV4-5NF (5’-AAGCTCGTAGTTGAATTTCG-3’) (Shi et al., 2019), and the following program: 3 min of denaturation at 95 °C, 30 cycles of 30 s at 95 °C, 30 s for annealing at 55 °C, and 45 s for elongation at 72 °C, and a final extension at 72 °C for 10 min. The resulted PCR products were extracted from a 2% agarose gel and further purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA) and quantified using QuantiFluor™-ST (Promega, Madison, WI, USA) according to the manufacturer’s protocol.

Illumina MiSeq DNA sequencing

Purified barcoded amplicons were pooled in equimolar concentrations and paired-end sequenced on an Illumina MiSeq PE300 platform/NovaSeq PE250 platform (Illumina, San Diego, CA, USA) according to the standard protocols by Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China). The raw reads were deposited into the NCBI Sequence Read Archive (SRA) database (Accession Number: PRJNA694003).

Processing of sequencing data

The raw sequencing reads were demultiplexed, quality-filtered by fastp version 0.20.0 (Chen et al., 2018) and merged by FLASH version 1.2.7(Magoč & Salzberg, 2011). Microbial community sequencing was conducted by Shanghai Majorbio Bio-pharm Technology using Illumina-MiSeq sequencing platform. The data were analyzed on a free online platform of Majorbio Cloud Platform (www.majorbio.com).

Operational taxonomic units (OTUs) were clustered with 97% similarity cutoff using UPARSE version 7.1, and chimeric sequences were identified and removed (Stackebrandt & Goebel, 1994; Edgar, 2013). The taxonomy of each OTU representative sequence was analyzed by RDP Classifier version 2.2 against the maarjam081/AM database using confidence threshold of 70% (Wang et al., 2007).

The raw sequence data were deposited at NCBI, in the SRA database with the following accession: PRJNA694003.

Data analysis

Total soil carbon, nitrogen, and sulphur were determined by an elemental analyser (GC IsolinkFlash 2000; Thermo Scientific, Waltham, MA, USA) analyzer. The concentration of total C, N, and S in soil were 6.96%, 0.55%, and 0.05%, respectively. Meanwhile, the C/N was 12.62. At the same time, we found that the soil temperature was increased 1.4 °C. The dynamic range of the soil temperature was increased from 0.6 °C to 2.4 °C in different elevations. The soil moisture decreased 0.07 m³. The soil moisture increased 0.15 m³ at 4,170 m and decreased 0.11, 0.09, and 0.25 m³ at 3,000, 3,500, and 3,800 m.

The community was expressed by AMF richness, relative abundance and occurrence frequency in different elevations. AMF richness was calculated by the number of OTUs. The relative abundance of AM fungal genus was calculated as the percentage of the sequence number of OTUs in each genus divided by the total sequence number of OTUs in all genera at this altitude. The occurrence frequency of AM fungal genus was defined as the percentage of the number of samples where this genus observed to the number of all samples in this genus. The rate of decrease = (the number of OTUs in CK − the number of OTUs in OTC)/the number of OTUs in CK * 100%. The rate of increase = (the number of OTUs in OTC − the number of OTUs in CK)/the number of OTUs in CK * 100%.

AMF alpha diversity in different elevations were expressed and plotted by the index of Shannon, Ace, and Simpson evenness at the level of OTU by SPASS, Excel and Origin, respectively. The differences of AMF richness and AMF diversity in different elevations were analyzed by two-way ANOVA analysis and Duncan in SPSS 19.0 (Shi et al., 2019). We analyze the impact of environmental factors on AMF community after warming through RDA analysis on the platform of the Majorbio Cloud Platform (www.majorbio.com). We standardize the data by flattening according to the minimum number of sample sequences on the platform of the Majorbio Cloud Platform (www.majorbio.com). The data of the percentage of relative abundance and occurrence frequency were subjected to square root transformation before analysing and comparing (Shi et al., 2007).

AMF richness at the level of OTU

Warming increased AMF richness at the level of OTU from 36 to 45.67 with the increase of 26.86%, among them, AMF richness of shared was 28 OTUs, which was 77.78% and 61.31% in CK and OTC, respectively (Fig. 1). In CK, there were 8 unique OTUs, which was 22.22% of the total in CK. AMF richness of shared was 3.5 times to CK solely. In OTC, there were 17.67 unique OTUs, which was 38.69% of the total in OTC. AMF richness of shared was 1.58 times to AMF richness in OTC solely. AMF richness increased but has no significant effects after warming by ANOVA analysis (P = 0.052).

AMF diversity indices at the level of OTU based on different elevations

There were dynamic influences of warming on AMF OTU richness with elevations. OTU richness displayed an upward trend at 3,000, 3,500 and 3,800 m but then decreased at 4,170 m after warming (Fig. 2A). The highest AMF richness occurred at 3,500 m and AMF OTU richness in OTC is greater than that in CK at the elevations of 3,000, 3,500, and 3,800 m, but it was opposite at 4,170 m. That was, AMF richness was lower at the higher altitude after warming. Moreover, elevation had extremely significant effects on AMF richness, which increased significantly at 3,500 m (P < 0.001). The interaction of elevations and warming also had a significant effect on AMF richness (P = 0.029). The Shannon index has the same tendency to AMF richness (Fig. 2B). At 3,000, 35,00, and 3,800 m, the Shannon index in OTC were higher than that in CK, but showed an opposite trend at 4,170 m and none of them are significant. The Simpson evenness index had a similar trend to the Shannon index at 3,000 m (Fig. 2C).

At 3,000 and 3,500 m, the Ace diversity index in OTC were higher than that in CK, but it was opposite at 3,800 and 4,170 m (Fig. 2D). Elevation had significant effects on the Ace index (P = 0.045).

The influence of warming on AMF community based on different elevations

Among the genera of Ambispora, Unclassified (Archaeosporales order), and Paraglomus, AMF richness of CK was identical with OTC (Table 2). The largest change in AMF richness was Glomus, which increased from 13.67 to 20.33 after warming. For Acaulospora, AMF richness was increased from 10.67 OTUs to 12.33 OTUs. The smallest change in AMF richness was Archaeospora, which increased from 4 to 4.33.

In addition, there was a downward trend at 4,170 m and the decline rate was 100% in Ambispora. However, there was an increasing trend at 3,500 and 3,800 m. For Archaeospora, AMF richness was increased at 3,800 m, but decreased at 3,500 and 4,170 m. For Acaulospora, AMF richness trended to increase 7.85% and 26.38% at 3,000 and 3,500 m, respectively. As for 4,170 m, it decreased 42.86%. For Glomus, AMF richness was increased 116.5%, 58.53%, and 38.83% at 3,000, 3,500, and 3,800 m, respectively, and then decreased at 4,170 m, AMF richness of Paraglomus increased at 3,000 m and decreased at 4,170 m. Moreover, the rate of increase on AMF richness at 3,000 m was the same as the rate of decrease at 4,170 m.

The beta-diversity of the AMF community was determined by the Bray–Curtis method (Fig. 3). In the 2-dimensional NMDS plots, soil samples collected from the four different elevations and the two different treatment separated from each other indicating a divergence of the warming treatment. To test the significance, an ANOSIM based on the Bray–Curtis distance showed dissimilarities of the AMF community at the OTU level among the four different elevations and two different treatments (P = 0.001).

The relative abundance and occurrence frequency of AMF

For four different elevations, the relative abundance of Acaulospora and Glomus were the largest before and after warming (Table 3). The relative abundance of Acaulospora and Glomus were opposite in four different elevations. At 3,000 m, the relative abundance of other genera showed an increasing trend after warming, except Glomus was decreased from 79% to 65%. At 3,500 m, the relative abundance of Ambispora and Glomus showed an increasing trend after warming, but Archaeospora and Acaulospora decreased. At 3,800 m, the relative abundance of Glomus showed an increasing trend but the relative abundance of Ambispora, Archaeospora, and Acaulospora showed a decreasing trend after warming. As for 4,170 m, the relative abundance of Archaeospora and Acaulospora increased but Ambispora, Glomus and Paraglomus decreased. The relative abundance of Glomus decreased, but Acaulospora increased after warming at 3,000 and 4,170 m.

For the four different elevations, Acaulospora was always present (Table 4), as was the occurrence of Glomus, except at 4,170 m in the OTC treatment. In CK, the occurrence frequency of Acaulospora was the same as that in OTC at different elevations, which seemed that warming had no affect on them. The occurrence frequency of Ambispora and Archaeospora varied at three elevations, Paraglomus varied at two elevations and Glomus varied at only one elevations. For different elevations, the occurrence frequency of Paraglomus showed a tendency of increasing from 0 to 33.33% at 3,000 m but opposite at 4,170 m. The occurrence frequency of Ambispora and Archaeospora decreased at 4,170 m. But at 3,500 m, the tendency of Ambispora and Archaeospora were opposite.

The influence of soil factors on AMF community by warming

For the four different elevations, RDA1 explained 59.16%, 69.55%, 95.26%, and 95.35% at 3,000, 3,500, 3,800, and 4,170 m of the community structure, respectively. RDA2 explained 0.91%, 2.34%, 1.52%, and 0.56% at 3,000, 3,500, 3,800, and 4,170 m, respectively (Figs. 4A–4D). RDA1 increased from 59.16% to 95.35% with the elevation increased. As the elevations increased, the influence of C, N, S, and C/N were different. C, N, and C/N were positively correlated to RDA1 and RDA2. Sulfur (S) was negatively correlated to RDA1 but positively correlated to RDA2 at 3,000 m (Fig. 4A). At 3,500 m, C, N, and S were all positively correlated to RDA2 but negatively correlated to RDA1 (Fig. 4B). C/N was negatively correlated to RDA2 but positively correlated to RDA1. C and N were negatively correlated to RDA1 and RDA2. S and C/N were negatively correlated to RDA2 but positively correlated to RDA1 at 3,800 m (Fig. 4C). At 4,170 m, C, N, and S were negatively correlated to RDA1 and RDA2. C/N was negatively correlated to RDA2 but positively correlated to RDA1 (Fig. 4D).