Rapid response of arbuscular mycorrhizal fungal communities to short-term fertilization in an alpine grassland on the Qinghai-Tibet Plateau

Site selection and soil sampling

The fertilization addition experiment started in 2011 (37°37′N, 101°12′E; 3,220 m), at the Haibei Alpine Grassland Ecosystem Experimental Station (QTP station), which was established in 1976 (northeast of QTP; 37°29′–37°45′N, 101°12′–101°23′E). The Haibei Alpine Grassland Ecosystem Experimental Station provided us with the necessary permits to carry out this work. The field sites were strictly vetted based on our selection standards: Mat-Cryic Cambisols soil type, vegetation type, and environmental conditions. The area has a typical plateau/continental climate, with a long winter. The average annual temperature is −1.7 °C, with a maximum temperature of 27.6 °C and a minimum temperature of −37.1 °C.

Soils were collected on the 12^th of August, 2014, at QTP Station. The experimental plots were constructed in randomized block design and the plot was a 6 × 6 m square. Our study included four treatments, with five replicates each: Control (without fertilization); N (100 kg N ha⁻¹ yr⁻¹); P (50 kg P ha⁻¹ yr⁻¹); and NP (100 kg N ha⁻¹ yr⁻¹ and 50 kg P ha⁻¹ yr⁻¹) (Table S1). Our fertilizer amendment amounts fall within an intermediate range of what other studies report (Sun et al., 2015; Zeng et al., 2015). The fertilizers were urea and Ca(H₂PO₄)₂, applied on the first days of June, July and August (i.e. during the main growing season). In each plot, soils were collected from four corners of the plot (1 m from edge of plot to avoid edge effects and at a depth of 0–10 cm) and then mixed together into one sample. The soils were kept in a cooler and shipped refrigerated to the lab as quickly as possible. The samples were completely mixed within each bag and sieved to remove stones, and then divided into three parts: one part was for biogeochemical analysis and was stored at 4 °C; the second part was stored at −20 °C for DNA extraction and the third part was placed in long-term storage at −40 °C.

Vegetation and soil properties analyses

Aboveground Net Primary Productivity (ANPP) were assessed in four 0.5 × 0.5 m areas, at the corners of soil plots. Aboveground portion of biomass was collected from a plot using clipper. Plants were grouped visually into four categories: Gramineae, Sedge, Legume and Forb. In each plot, the roots from three replicate soil cores which were collected using a 7 cm diameter auger were used to estimate belowground root biomass. The dry biomass for aboveground sections and belowground roots (after washing) were calculated after 48 h in the 65 °C oven. Soil pH was measured after shaking a soil water suspension (1:5 wt/vol) for 30 min and soil moisture (SM) was measured gravimetrically. The classical methods were applied for measuring soil available phosphorus (AP, Ståhlberg, 1980), total phosphorus (TP, Bowman, 1988), total carbon (TC) and total nitrogen (TN) (Walkley & Black, 1934). Soil dissolved organic C (DOC) and total dissolved N (TDN) and mineral nitrogen were extracted by adding 50 ml of 0.5 M K₂SO₄ to 10 g fresh soil, shaking for 1 h and then vacuum filtering through glass fiber filters (Fisher G4, 1.2 μm pore space). Ammonium (NH₄⁺) and nitrate (NO₃⁻) contents in the extracts were determined colorimetrically by automated segmented flow analysis (Bran + Luebbe AAIII, Germany) using the salicylate/dichloroiso- cyanuric acid and cadmium column/sulphanilamide reduction methods, respectively. DOC and TDN were determined using a TOC-TN analyzer (Shimadzu, Kyoto, Japan). Dissolved organic N (DON) was calculated as follows: DON = TDN − NH₄⁺ − NO₃⁻. Biogeochemical data are shown in Table S2.

Soil DNA extraction

DNA extractions were carried out on 0.5 g soil according to the manufacturer’s instructions (FastDNA® SPIN Kit for soil, MP Biomedicals, Santa Ana, CA). The extracted DNA was dissolved in 50 μl TE buffer, quantified by NanoDrop ND-1000 (Thermo Scientific, USA) and stored at −20 °C.

PCR and amplicon library preparation

Primers AMV4.5NF and AMDGR (Lumini et al., 2010) were used to amplify soil 18S rRNA gene fragments for the Illumina MiSeq platform (PE 300) at Novogene (Beijing, China). PCR was carried out in 50-μl reaction mixtures containing each deoxynucleoside triphosphate at a concentration of 1.25 mM, 1 μl of forward and reverse primers (20 μM), 2 U of Taq DNA polymerase (TaKaRa, Japan), and 50 ng of DNA. The following cycling parameters were used: 35 cycles (95 °C for 45 s, 56 °C for 45 s, and 72 °C for 45 s) were performed with a final extension at 72 °C for 10 min. Triplicate reaction mixtures per sample were pooled together and purified using an agarose gel DNA purification kit (TaKaRa), and quantified using NanoDrop ND-1000 (Thermo Scientific, Waltam, MA, USA) ranging from 44.2 to 82.7 ng/μl. The bar-coded PCR products were pooled in equimolar amounts (10 pg for each sample) before sequencing.

Processing of sequence data

Sequences were merged by FLASH (Magoc & Salzberg, 2011) and then processed using Quantitative Insights Into Microbial Ecology (QIIME; http://www.qiime.org/) (Caporaso et al., 2010). Poor-quality sequences (below an average quality score of 25) and short sequences (< 200 bp) were removed. Sequences were clustered into Operational Taxonomic Units (OTUs) using a 97% identity threshold (default QIIME settings) by UCLUST (Edgar, 2010) and all singleton OTUs were deleted. Chimera filtering was also performed to remove sequencing errors with USEARCH tool (version 1.8.0) in QIIME. The most abundant sequence within each cluster was selected as the representative sequence for that OTU. The representative sequences were checked against the MaarjAM AMF database (Öpik et al., 2010; http://maarjam.botany.ut.ee/). We rarified the abundance matrix to 1,000 sequences per sample to obtain normalized relative abundances.

Statistical analysis

Phylogenetic diversities (PD) were estimated by Faith’s index (Faith, 1992), providing an integrated index of the phylogenetic breadth across taxonomic levels. Pearson correlations were calculated between AMF alpha-diversity and biogeochemical properties using SPSS 20 for Windows. A neighbor-joining tree (MEGA 6; Tamura et al., 2013) was built for phylogenetic assignment by aligning dominant AMF OTU sequences (relative abundance > 1%) with known reference sequences. The response ratio (RR) was used to quantify significant responses of OTUs to fertilization. Only those OTUs detected in more than 3 replicates of each treatment were analyzed. The $\text{95}\%\text{\hspace{0.17em}confidence\hspace{0.17em}interval\hspace{0.17em}}\left(\text{CI}\right)=r{r}_i\pm 1.96\times \sqrt{V_i}$, where $r{r}_i=\text{ln}\left({\overline{x}}_i/{\overline{y}}_i\right)$ (i = 1 … n), $\overline{x}$ is the mean OTU read number in fertilization samples and $\overline{y}$ is the mean OTU read number in control samples; the variance (V_i) is, $V_i=\frac{s_{x_i}{}^2}{m_{x_i}{\overline{x}}_i{}^2}+\frac{s_{y_i}{}^2}{m_{y_i}{\overline{y}}_i{}^2}$ (i = 1 … n), where s is the SD of OTU i in control and fertilization samples and m is the abundance of OTU i in control and fertilization samples. The rr-significance value was calculated using the equation: $rr\text{-}significance=\left(r{r}_i+1.96\sqrt{V_i}\right)\times \left(r{r}_i-1.96\sqrt{V_i}\right)$. If the 95% CI of a response variable overlaps with zero (rr-significance < 0), the RR at fertilization is not significantly different from the control. If rr-significance > 0, the RR is significant (Luo, Hui & Zhang, 2006; Xiang et al., 2015). Non-metric multidimensional scaling (NMDS) based on Bray-Curtis dissimilarity (calculated from the relative abundance matrix), and Analyses of Similarities (ANOSIM) (Clarke, 1993) were performed to compare community composition in different treatments in the vegan package 2.0–10 (Dixon, 2003) of R v.3.1.0 (R Development Core Team, Vienna, Austria). Mantel tests were used to identify biogeochemical factors that significantly correlated with community composition. Multivariate regression trees (MRT) were constructed to display important relationships between AMF community and biogeochemical variables in the mvpart package 1.2–6 (De’ath, 2002).

General sequencing information

We obtained a total of 380,203 quality-filtered fungal sequences across all soil samples for the primer pair AMV4.5NF/AMDGR: 12.14% were classified as Glomeromycota (AMF). A total of 46,140 quality-filtered AMF sequences were identified across all the entire data set, ranging from 1,088–4,313 sequence reads per sample (mean = 2,307). The neighbor-joining tree showed that the dominant OTUs (relative abundance > 1%) grouped into Claroideoglomeraceae (30.1%), Gigasporaceae (23.9%), Glomeraceae (23.4%), Diversisporaceae (12.6%) and Acaulosporaceae (9.1%) families (Figs. S1 and S2). In addition, Ambisporaceae, Archaeosporaceae, Pacisporaceae were present in most soils but at relatively low abundances, indicating a good coverage of the Glomeromycota.

Arbuscular mycorrhizal fungal alpha-diversity

After rarefaction to an equal sequencing depth per sample (1,000), we found a total of 231 distinct AMF OTUs across all samples (97% similarity), 35.5% of which (82) were found in all treatments. There were 132, 131, 143 and 169 AMF OTUs in Control, N, P and NP treatments, respectively. In addition, there were unique OTUs within each treatment, especially for NP plots, where the unique OTUs accounted for 19.5% of the 169 distinct OTUs observed (Fig. 1). Soil AMF alpha-diversity (i.e. OTU richness and phylogenetic diversity) was calculated at a depth of 1,000 randomly selected sequences per sample. Compared to control soils, NP addition significantly increased AMF OTU richness and phylogenetic diversity in this study (Fig. 2). Of all the biogeochemical characteristics examined, both OTU richness and phylogenetic diversity were positively correlated with NO₃⁻, DON, AP, graminoid biomass, forb biomass and TB and negatively correlated with legume biomass (P < 0.05 in all cases; Table S3).

Arbuscular mycorrhizal fungal community composition

Compared to control plots, NP addition significantly decreased the relative abundance of Gigasporaceae and increased the relative abundance of Glomeraceae. The relative abundance of Diversisporaceae showed significant increase and the relative abundance of Acaulosporaceae showed significant decrease in N, P and NP plots relative to control soils (Fig. 3). RR were calculated in order to identify significant differences across control and fertilization samples (Fig. 4). In the family Diversisporaceae, all OTUs showed significant increases in abundance, and for Acaulosporaceae, all the OTUs showed significant decreases in abundance in fertilization plots. In the family Gigasporaceae, all OTUs in P and NP plots showed negative responses relative to the control, and for the family Glomeraceae, all OTUs showed positive responses in NP plots relative to the control.

AMF community composition was significantly influenced by short-term fertilization (Fig. 5; Table S4). Mantel tests showed that the AMF community composition was significantly correlated with AP, NO₃⁻, TP, graminoid biomass, legume biomass, forb biomass and total aboveground biomass (TB) (P < 0.05 in all cases; Table S5). In addition, MRT analysis was used to investigate the effects of biogeochemical variables on AMF communities. The MRT model explained 37.0% of detected variation in AMF community composition. The important factors were AP, TB and legume biomass, which accounted for 20.5, 8.4 and 7.1% variation of AMF community composition, respectively (Table S6). Both AP and TB showed significant Pearson correlations with the relative abundance of Claroideoglomeraceae, Gigasporaceae and Acaulosporaceae; Both soil DON and forb biomass showed significant Pearson correlations with the relative abundance of Glomeraceae and Acaulosporaceae; TP showed significant Pearson correlations with the relative abundance of Gigasporaceae and legume biomass showed significant Pearson correlations with the relative abundance of Glomeraceae (P < 0.05 in all cases; Table S7).