The Hailuogou glacier (29°24′N, 101°59′E) is located in Mount Gongga, which is the highest peak in the Hengduan Mountain region in the eastern part of and the south-eastern edge of the Tibetan Plateau. As a typical monsoonal temperate glacier, the Hailuogou glacier is the largest in the basin and has a total area of c. 25 km² with c. 13 km long (Li et al., 2010). Primary vegetation succession has developed on this chronosequence (Li et al., 2010). Mean annual temperature is 4.2°C. The diurnal temperature range is 20.4°C. Mean annual precipitation is 1,947 mm.

Soil samples were collected from six sites along the chronosequence, ranging from <100 to 2,000 m from the glacier terminus (Figure 1). Correspondingly, these sites have deglaciated for c. 2, 18, 38, 68, 110, and 188 years, resulting in vegetation that covers 1, 4, 70, 75, 85, and 90% of the surface, respectively. A detailed description of the sampling sites is listed in Table 1 and Table S1. At each site the samples were collected from six replicate plots (10 × 10 m). Each soil sample was immediately placed into a sterile sample bag and sealed, and transported to the lab in an ice-box. One portion was stored at −80°C for DNA extraction while the other half were stored at 4°C for the soil analysis.

The pH-value was determined at a soil to water ratio of 1:5 (wt/vol). Total nitrogen (TN) was determined using a modified Kjeldahl method. Ammonium (NH4+–N) and nitrate (NO3−–N) were measured using a continuous flow analyzer (Skalar SA 1000, Breda, The Netherlands). The soil samples for analysis of the dissolved organic carbon (DOC) were prepared by mixture with 0.5 M K₂SO₄ in a centrifuge tube, then shaking for 1 h on a reciprocal shaker and centrifuging at 13,000 rpm for 30 min at 4°C. The supernatant was filtered through a 0.45-μm glass fiber filter. All the extracts were analyzed by a continuous flow analyzer of San++ (Skalar SA 1000, Breda, The Netherlands).

Soil microbial biomass carbon (SMBC) and soil microbial biomass nitrogen (SMBN) were determined by the chloroform fumigation extraction method after 14 days of conditioning at 50% of their total water holding capacity under 25°C. This was followed by a 0.5 M K₂SO₄ extraction method for both the non-fumigated and fumigated samples (Brookes et al., 1985; Vance et al., 1987). The SMBC extraction was analyzed by a TOC-V WP (SHIMADZU, Japan) and the SMBN extraction was digested and then analyzed by a San++ continuous flow analyzer (Skalar SA 1000, Breda, The Netherlands).

Total DNA for the metabarcoding analyses was extracted from c. 0.5 g of frozen soil with the UltraClean Soil DNA extraction kit (MoBio, Carlsbad, USA), according to the manufacturer's instructions. The obtained DNA was quantified with a NanoDrop device and stored at −20°C. An aliquot of 50-ng purified DNA from each sample was used as template for amplification. A set of primers was designed by adding a Roche 454 “A” pyrosequencing adapter, and a unique 8-bp barcode sequence to the forward primer ITS1F (5′-CTTGGTCATTTAGAGGAAGTAA-3′) and the reverse primer ITS2 (5′-GCTGCGTTCTTCATCGATGC-3′) (White et al., 1990) to generate PCR ITS1 region fragments of c. 400 bp. Each sample was amplified in triplicate with 50 μl reactions under the following conditions: 94°C for 4 min, 30 cycles of 30 s at 94°C (denaturation), 50°C for 1 min (annealing), and 72°C for 90 s (extension), followed by 10 min at 72°C. PCR products were pooled together and purified by the Agarose Gel DNA purification kit (TaKaRa, China) and quantified with the NanoDrop device. Samples were run on a ROCHE 454 FLX+ platform (Roche, Basel, Switzerland) at the National Human Genome Center of China at Shanghai, China, according to the manufacturer's instruction manual. These raw data are deposited in the NCBI Sequence Read Archive with the accession number SRP093928.

The raw sequence data were de-multiplexed using the QIIME toolkit (Caporaso et al., 2010). Reads with a quality score of less than 20 were filtered out (split_libraries_fastq.py; QIIME) (quality score <20, homopolymer runs of >6 nt and length <200 nt), and the remaining reads were assigned to samples according to their corresponding barcode. The extracted sequences were binned into operational taxonomic units (OTUs) at 97% identity using UPARSE (Edgar, 2013); species-level similarity was denoted by a closed-reference OTU that used the UNITE ITS database (Abarenkov et al., 2010) as a reference (pick_otus.py; QIIME). One representative of each OTU was selected for the downstream analyses, and taxonomy was assigned to each representative (pick_rep_set.py, assign_taxonomy.py; QIIME). The α-diversity of each sample was calculated using the Shannon metric (alpha_diversity.py; QIIME).

All singleton OTUs from the downstream analyses were excluded because whether singleton data represents true diversity or sequencing errors is an unsettled issue (Tedersoo et al., 2010). Evenness and Shannon diversity indices (H′) were estimated for total fungi as well as their most abundant taxa based on OTU abundance matrices rarefied to the lowest sequence numbers. Chao and ACE were calculated at the 97% threshold. Sampling effort was estimated using Good's coverage.

To evaluate the effects of distance from the glacier terminus on OTU richness, Shannon diversity, and evenness, we fitted linear and cubic regressions. In these analyses we treated distance from the glacier terminus as a continuous (linear, cubic regression). The effect of distance on the relative abundances of the most abundant and frequency taxa were assessed by linear and quadratic regressions. All statistical analyses were performed in R.

Principal coordinate analyses (PCoA) were used to compare dissimilarities among the fungal communities across succession distance, as based on the Bray-Curtis dissimilarity. To determine the relationship between succession distance and Bray-Curtis distances, scores for the first two principal coordinates against succession distance were regressed. Mantel test is multivariate correlation analysis by means of comparison between dissimilarity matrices of same dimension. Mantel tests were performed to examine community turnover along environmental and succession gradients. The overall effect of glacier retreating on fungal community was determined using and permutational MANOVAs (PERMANOVAs). These above analyses were performed in the vegan package in R.

β-diversity is useful to understand patterns of species diversity across various spatial scales and provides critical insights into the role of deterministic and stochastic processes in shaping community compositions and structure (Chase, 2007, 2010; Chase and Myers, 2011; Zhou et al., 2014). Low dissimilarity among communities that are otherwise identical in environmental conditions would imply a predominant role for deterministic assembly, whereas high dissimilarity would suggest a large role of stochastic assembly (Chase, 2007, 2010). However, the measures of β-diversity are dependent on both α- and γ-diversity (Chase et al., 2011). It is not clear whether a change in β-diversity is due to the differences in the underlying assembly processes that generate β-diversity or to the differences in α- and γ-diversity (Chase et al., 2011). β-Diversity null deviation approach uses a null model to create stochastically assembled communities from the regional species pool to determine the degree to which the observed β-diversity patterns deviate from stochastic assembly(Chase and Myers, 2011; Tucker et al., 2016). Hence, this approach can assess changes in β-diversity that result from the relative influences of deterministic and stochastic processes, and not from changes in α-diversity. We measured the null deviation as the relative difference of the observed β-diversity from the null-model β-diversity, (β_obs−β_null)/β_null, where β-diversity was measured as the Sorenson-Czekanowski dissimilarity. For each sample, the expected β-diversity under the null model was calculated from 999 stochastically assembled communities. The γ-diversity was calculated from the total number of species detected in the “species pool” from all sites. To test for treatment differences in the null deviation, we conducted permutation tests by first randomly permuting treatment labels, then re-simulating the null models and re-calculating the null deviations for each of the 999 permutations.

Standard effect size (SES), another similar quantitative index, was introduced into this study to measure the influence of deterministic factors on community composition and abundance over time. β-Deviation was calculated as the observed β-diversity subtracted the mean of the null distribution of β-diversity values, then divided by the standard deviation of this distribution (Kraft et al., 2011).

After quality control (QC), 52% of fungal ITS-1 sequences were removed through each step of our pipeline and only 145,450 sequences from 299,945 were remained (Table S2). Sequence loss through each step of our pipeline is summarized in Table S2. Clustering of the sequences resulted in a total of 611 OTUs. The coverage of each sample ranged from 0.92 to 0.99 (Table 1), suggested that these libraries detected a large majority of the fungal diversity in the samples used in our study. Soil fungal communities were strongly dominated by the diverse Ascomycota and Basidiomycota (47.0 and 48.3% of total sequences, respectively). Surprisingly, a relatively higher component of Rozellomycota (Cryptomycota) was detected for 4.49% of total sequences and 38 OTUs (6.2%), compared with that of Chytridiomycota and Glomeromycota, which was constituted of 0.2 and 0.01% of the total sequences, respectively. The sequencing data suggested that the fungal community was uneven and that 10 dominant OTUs (abundance >1%) out of the 611 OTUs accounted for 70% of the total reads. Among those 10 OTUs, three most abundance OTUs were assigned to Cysrodendron sp. (25.5% of total sequences), and followed by Cryptococcus terricola (16.4%), and Scleroderma areolatum (15.9%).

A positive linear relationship was found between the relative abundance of Rozellomycota and nitrate concentration (Linear regression, R² = 0.34, P < 0.001). No significant relationship was detected between nitrate concentration and Ascomycota and Basidiomycota.

The fungal indices for Shannon diversity and evenness showed unimodal patterns over distance from the glacier terminus (Figures 2A,B), while Chao1 richness and ACE showed no significant response trend (Figures 2C,D). Shannon diversity and evenness peaked at a distance of 370 m and declined afterwards. The Shannon diversity index of Basidiomycota significantly decreased with distance (Figure 3A, linear regression, P < 0.05), while the pattern of Rozellomycota diversity response to succession distance was unimodal (Figure 3B, quadratic regression, P < 0.05). No significant response pattern was detected between Ascomycota diversity and succession distance (Figure 3C, linear regression and quadratic regression, P > 0.05).

To verify the differences observed in the fungal communities along chronosequence, the relative abundances of the different classes were compared. The relative abundance of Rozellomycota (Cryptomycota) significantly increased in younger developing soil but declined after 900 m (Figure 4A). The relative proportion of Tremellomycetes increased (Figure 4B), whereas those of Sordariomycetes (Figure 4C), and Leotiomycetes decreased with distance (Figure 4D). Agricomycetes did not significant response to distance (Figure 4E).

Ecosystem processes and successional dynamics are probably driven by those community members that occur most frequently. There were eight OTUs that occurred at all the six sites (Figure 5). Interestingly, the relative abundance of OTU2 (C. terricola) significantly increased with increased distance from the glacier terminus (Linear regression, P < 0.001; Figure S1), but its relative abundance was negatively related to fungal diversity (Figure S2). However, the relative abundance of OTU65 (Tolypocladium tundrense) decreased with distance from the glacier terminus (Linear regression, P < 0.001).

Principal coordinate analyses (PCoA) and PERMANOVA analyses revealed significant differences in the fungal communities along the successional gradient (Figure 6, PERMANOVA, pseudo-F = 3.82, P < 0.001). The fungal communities from younger soils (distance: <100 and 200 m) were significantly differentiated to samples from distance of 370 m, indicated the soils are developing. Whereas, the communities from the older soils (distance: 900–2,000 m) were more similar to each other, indicated the soils are developed (Figure 6). The scores for axis 1 representing 21.6% of the variability increased substantially along with distance from the glacier terminus (t = 2.43, P < 0.05). Treating distance from the glacier as a categorical variable, the axis 1 scores differed among distance classes [ANOVA: F_{(1, 32)} = 5.89, P < 0.05]. This was mainly attributable to the axis 1 scores between the younger soils and older soils (Tukey's HSD pairwise comparison at an alpha level = 0.05). The scores for axis 2 representing 10.5% of the variability also tended to increase with distance from the glacier terminus (t = 2.59, P = 0.02).

Mantel test results showed that distance (r = 0.26, P = 0.003), altitude (r = 0.19, P = 0.02), SMBN (r = 0.37, P = 0.016), SMBC (r = 0.34, P = 0.016), and NH4+–N (r = 0.34, P = 0.06) correlated with fungal community.

The null deviation approach can provide powerful insights into community assembly mechanisms (Chase and Myers, 2011). A null deviation close to zero suggests that stochastic processes are more important in shaping community structure, whereas higher positive or negative null deviations suggest that deterministic processes are more important (Chase, 2007, 2010; Ferrenberg et al., 2013). The relative contribution of both processes was dynamic along the successional stages tested (Figure 7). In the initial stage, the fungal communities significantly deviated from the stochastic assembly model (relative null deviation = 0.67–0.74) more than in later stages (relative null deviation = 0.45–0.52) (P < 0.05). Importantly, the intermediate stages and later stages showed a relatively consistent deviation from the stochastic assembly model, with no significant changes in the null deviation value between samples (Figure 7). Also, SES, which is used to measure the influence of deterministic factors on community composition and abundance (Kraft et al., 2011), was higher for the communities at the early stages than intermediate and later stages (Figure S3).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.