The LHG (39°26.4′N, 96°32.5′E), with a length of 10.1 km and an area of 21.9 km², is the largest glacier in the Qilian Mountains, in the northeastern Tibetan Plateau (Sun et al., 2012). It is a typical valley glacier, with elevations ranging from 4,200 m a.s.l at its terminus to 5,481 m a.s.l at its top. The region has a typical continental climate affected by desertification in central Asia, with an average temperature above 0°C in summer and long periods below freezing in winter. Precipitation is mainly concentrated in May-September, primarily affected by westerly winds (Wei et al., 2017).

In August 2016, we randomly selected four superficially deposited cryoconite mounds, which piled up due to glacial melting within the ablation zone for collecting samples (Supplementary Figure 1). Four sampling sites were 100–400 m apart, covering a large area in order to represent the glacier’s features. To ensure that collecting was made at the same locations, we determined the coordinates of each site using GPS and mounted bamboo sticks in the ice. Samples were collected from each of the observation sites on the 1st, 4th, 8th, 26th, and 31st of August, which are referred as day 1, 4, 8, 26, and 31 hereafter). The large gap between the third and fourth sampling days was due to heavy snowfall and flooding that occurred after 10th August. Before sampling, we scraped off the surface layer (top 1 cm) of cryoconite. Approximately 200 g of cryoconite debris were collected using a pre-cleaned stainless-steel spoon with polyethylene gloves worn all the time. The sample was placed into sterile 250 ml Nalgene ^® plastic bottles that were first cleaned with tap water and then rinsed three times with Milli-Q water, followed by sterilization for 15 min at 120°C in steam autoclaving and drying at 55°C. Samples were transported in a dark container within 4 h to the Qilian Shan Station of Glaciology and Ecologic Environment, Chinese Academy of Sciences. All samples were kept frozen at −20°C until laboratory analysis.

The total organic carbon (TOC) contents of the samples were determined using a TOC-L (Shimadzu Corp., Kyoto, Japan) on lyophilized samples. Total nitrogen (TN) and total carbon (TC) were measured using an elemental analyzer (Vario MAX, Elementar, Germany). The elemental concentrations are reported as the percentage of dry weight (100°C for 24 h) of cryoconite material.

Cryoconite debris (0.5 g) was soaked in 1.5 mL of lysis buffer (20 mg mL^–1 Proteinase K, 0.1 M EDTA, and 10% SDS) and incubated at 55°C for 2 h before DNA extraction. DNA extraction was performed using a Fast DNA ^®SPIN Kit for Soil (MP Biomedicals, Santa Ana, CA, United States), according to the manufacturer’s instructions. The extracted DNA was eluted to 100 μL TE buffer. The raw DNA was purified using a 1% (w/v) agarose gel. DNA bands were removed from the gel and extracted using an Agarose Gel DNA purification kit (TaKaRa, Japan). The concentration of DNA extracts was measured by the NanoDrop 1000 spectrophotometer (Thermo-Scientific, Wilmington, DE, United States). The extracted DNA was stored at -80°C until amplification.

The V4 region of the bacterial 16S rRNA genes was amplified in triplicate by primer pair 515F (5′-GTGCCAGC MGCCGCGGTAA-3′) and 806R (5′- GGACTACHVGGGTW TCTAAT-3′) and sequenced on an Illumina MiSeq system (Illumina, Inc., San Diego, CA, United States) following the procedure by Caporaso et al. (2012). Briefly, PCR was performed under the condition of 94°C for 5 min, 30 cycles of 94°C for 30 s, 52°C for 30 s, 72°C for 30 s; followed by a final cycle of 10 min at 72°C. Each PCR reaction contained 25 μL 2x Premix Taq DNA polymerase (Takara Biotechnology, Dalian Co., Ltd., China), 3 μL DNA template (20 ng μL^–1), 1 μL each primer (10 μM), and 20 μL nuclease-free water. PCR products of the triplicate reaction samples were pooled and quantified by PicoGreen staining. The sample libraries were generated with purified PCR products. The MiSeq 500 cycles kit was used for 2 × 250 bp paired-end sequencing.

The “fastq_filter” command in USEARCH was used to trim and filter the raw sequence reads (Edgar, 2016). The trimmed and filtered reads were used to call amplicon sequence variants (ASVs) using the “-cluster_ASVs” and “-unoise3” commands in USEARCH, respectively (Edgar, 2016). Singleton ASVs with fewer than nine reads were removed using the default “-minsize” values. The “-usearch_global” command with a minimum sequence identity of 99% was used to map ASVs sequences. Chimeric ASVs were removed using the “uchime_ref” command in USEARCH. Representative sequences were taxonomically classified using the “sintax” command along with the RDP training set (version 16) with confidence scores above the 0.5 cutoff (Chiellini et al., 2019). Normalization across samples was performed using a variance stabilizing transformation with the DESeq2 package in R (McMurdie and Holmes, 2014). After the taxonomy assignment, ASVs associated with chloroplasts, archaea, and unclassified sequences were excluded from subsequent analyses.

Rare taxa were defined as the ASVs with mean relative abundance < 0.1% in all samples in one group, whereas abundant ASVs were defined as the ASVs with a relative abundance > 1% in one or more samples in one group (Ji et al., 2020; Zhang et al., 2020; Dong et al., 2021). Taxa that did not fall into either the abundant or rare categories were defined as the oscillating taxa.

Wilcoxon rank-sum test and Kruskal-Wallis test were used to determine whether there were any statistically significant differences among groups of environmental variables. The alpha-diversity indices (Shannon-Wiener and Pielou’s Evenness) were calculated using the “vegan” R package (Oksanen et al., 2010). Wilcoxon rank-sum test and Kruskal-Wallis test were used to determine whether there were any statistically significant differences among groups of alpha-diversity. The significance of dissimilarity of community composition among groups was tested using permutational multivariate analysis of variance (PERMANOVA) based on Bray-Curtis distance metrics with the “adonis” function in the R package “vegan” (Oksanen et al., 2010). Test results with P < 0.05 were considered as statistically significant. The correlations between bacterial diversity and environmental variables were calculated using Spearman’s rank correlation, as implemented in the “ppcor” R package. Mantel test based on Spearman’s rank correlations was performed to show the correlation between variations of microbial communities and environmental variables using the “vegan” package of R (Oksanen et al., 2010).

Co-occurrence networks were explored using Molecular Ecological Network Analysis Pipeline (MENA2)¹ based on the Random Matrix Theory (RMT) and Spearman correlation (Deng et al., 2012). The resulting correlation matrix was analyzed with the RMT-based network approach to determine the correlation threshold for network construction, and the same threshold was used for all networks, so the topological properties of all networks are comparable. Then global network properties were calculated including individual nodes’ centrality, degree, betweenness, and clustering coefficient. Integrating all the results obtained from MENA, co-occurrence network was visualized by Gephi.

The normalized stochasticity ratio (NST) based on the Bray–Curtis dissimilarity was calculated using the “NST” package in R to estimate the determinacy and stochasticity of the bacterial assembly processes with high accuracy and precision (Ning et al., 2019). The NST index used 50% as the boundary point between more deterministic (< 50%) and more stochastic (> 50%) assembly processes.

The TOC concentration ranged from 1.04 to 3.39%. The TC concentration ranged from 2.10 to 3.85%. The concentration of TOC and TC of samples after snowfalls were similar between the samples before and after snowfalls (Wilcoxon rank-sum test, P = 0.13, and P = 0.73, respectively). The TN concentrations ranged from 0.16 to 0.24%, and that before the snowfall was significantly higher than that after snowfall (Wilcoxon rank-sum test, P = 0.003). The C/N ratio ranged from 9.47 to 22.62, and the C/N ratio of samples after the snowfall was significantly higher than those before the snowfall (Wilcoxon rank-sum test, P = 0.02; Figure 1). All environmental characteristics showed no significant difference between sampling sites (Kruskal-Wallis test, all P > 0.05; Supplementary Figure 2).

After quality filtering and the removal of low-quality sequences, a total of 867,658 reads were obtained, and were clustered into 872 ASVs (99% sequence similarity). The most abundant bacterial phyla were Bacteroidetes, Cyanobacteria, Betaproteobacteria, Gammaproteobacteria, Actinobacteria, Alphaproteobacteria, Chloroflexi, and Firmicutes, which account for 86.2 and 95.0% of sequence abundance before and after snowfalls, respectively (Figure 2). The composition of the abundant sub-community was similar to the whole community. However, the rare and oscillating sub-communities differed from the abundant sub-communities. The rare and oscillating taxa contained more taxonomic groups than the abundant taxa. Bacteroidetes, Alphaproteobacteria, Verrucomicobia, Deltaproteobacteria, and Acidobacteria were more prevalent in rare and oscillating sub-communities, while Cyanobacteria and Betaproteobacteria were more prevalent in abundant sub-community (Figure 2).

In the whole community, the relative abundance of Betaproteobacteria, Chloroflexi, Alphaproteobacteria, Gemmatimonadetes, and Deltaproteobacteria were significantly higher before snowfall than after snowfalls (Wilcoxon rank-sum test, all P < 0.05; Supplementary Figure 3). In the abundant sub-community, the relative abundance of Bacteroidetes, Chloroflexi, Alphaproteobacteria, and Deltaproteobacteria were significantly higher before snowfall than after snowfalls (Wilcoxon rank-sum test, all P < 0.05), while Actinobacteria were significantly higher after snowfall than before snowfalls (Wilcoxon rank-sum test, P < 0.01). In the rare sub-community, the relative abundance of Alphaproteobacteria and Deltaproteobacteria were significantly higher before snowfall than after snowfalls (Wilcoxon rank-sum test, all P < 0.05). In the oscillating sub-community, the relative abundance of Gammaproteobacteria, Firmicutes, Gemmatimonadetes, and Deltaproteobacteria were significantly higher before snowfall than after snowfalls (Wilcoxon rank-sum test, all P < 0.05; Supplementary Figure 3).

The abundant taxa (34-48 ASVs) accounted for 62.1–94.5% of all sequences. The rare taxa (167-495 ASVs) only accounted for 2.6–13.6% of all sequences. The oscillating taxa (40-128 ASVs) accounted for 2.2–29.4% of all sequences (Supplementary Table 1). After snowfall disturbance, the proportion of abundant bacterial taxa did not change significantly, while the percentage of rare taxa significantly increased, and the percentage of oscillating taxa significantly decreased (Supplementary Figure 4). Approximately 59% of abundant bacterial taxa before snowfall were still abundant after the snowfall, and 8% became rare bacterial taxa, and 35% became oscillating bacterial taxa. After the snowfall disturbance, a large proportion (93%) of rare bacterial taxa still being the rare bacterial taxa, while 1% became abundant bacterial taxa, and 1% became oscillating bacterial taxa, and an additional 5% disappeared. Only 29% of oscillating bacterial taxa in samples before snowfall were still oscillating in bacterial taxa after the snowfall, 2% became abundant taxa, and 69% became rare bacterial taxa (Supplementary Table 2). The abundant bacterial taxa that came from the rare and oscillating bacterial taxa after snowfall accounted for 23 and 19% of the community. They are primarily taxonomically classified as Flavobacterium (Bacteroidetes), Massilia (Betaproteobacteria), Polaromonas (Betaproteobacteria), Psychrobacter (Gammaproteobacteria), Cellulomonas (Actinobacteria), and Carnobacterium (Firmicutes) (Supplementary Table 3).

The Shannon-Wiener diversity of the whole, abundant, rare, and oscillating community were 1.74–4.90, 1.47–3.67, 3.79–5.38, and 2.72–4.56, respectively. The Pielou’s Evenness of the whole, abundant, rare, and oscillating community were 0.02–0.22, 0.07–0.60, 0.46–0.68, and 0.33–0.85, respectively. The Shannon diversity indices of the rare sub-community were significantly higher than those of the abundant and oscillating sub-communities (Wilcoxon rank-sum test, all P < 0.001; Supplementary Figure 5). After the snowfall, the Shannon-Wiener and Pielou’s Evenness indices were significantly reduced in all subcommunities (Wilcoxon rank-sum test, all P < 0.05; Figure 3). The Shannon and Evenness indices positively correlated with the TN concentration (Spearman’s correlation: r = 0.45, P < 0.05; r = 0.45, P < 0.05; Table 1).

The PCoA plot showed clear segregation between the samples before and after snowfall (Figure 4). The result was also confirmed by the PERMANOVA analysis (PERMANOVA: F = 7.57, P < 0.001). Mantel tests revealed significant correlations of the community composition with the TN concentration (Mantel test: r = 0.45, P = 0.002; Table 2).

The co-occurrence network of bacterial communities was performed at the ASV level based on Spearman’s correlation relationships to examine the interactions between microorganisms. The co-occurrence network of bacterial communities was more complex after snowfall (Figure 5). Before the snowfall, the network contained 379 nodes which were connected by 690 edges, while after the snowfall the network contained 454 nodes and 726 edges. In addition, the networks before and after snowfall had 675 and 692 positive edges, respectively, but only 15 and 34 negative edges, respectively. The network after snowfall showed higher values of module number (68) and modularity (0.84) than the network before snowfall (31 and 0.82, respectively), indicating increased degree of modularity (Table 3). The nodes in the network before snowfall were composed of abundant ASVs (12.0%), rare ASVs (57.6%), and oscillating ASVs (30.4%). The nodes in the post snowfall network were composed of a lower proportion of abundant and oscillating ASVs (7.0 and 11.7%), but a higher proportion of rare ASVs (81.3%).

In the community before snowfall, the node degree and betweenness of the abundant bacterial were significantly higher than that of rare bacterial and oscillating taxa (Figure 6), indicating that more information may be passed through abundant bacterial taxa. While in the community post snowfall, no significant differences were detected between abundant, rare and, oscillating taxa.

NST was used to examine the relative contributions of stochasticity and determinism in shaping bacterial communities. The average NST values were 62 and 33% in the communities before and after snowfall (Figure 7). After the snowfall, the stochasticity processes significantly decreased for the whole community (Wilcoxon test, P < 0.001), suggesting the deterministic processes increased after snowfall disturbance.

Bacterial diversity reduced significantly after snowfall disturbance, indicating that the cryoconite microbial community may be subjected to greater environmental filtering after snowfall. This is mostly likely due to nitrogen limitation, as nitrogen was the only measured environmental factors whose concentration significantly reduced post snowfall. This was consistent with the positive correlative relationships between nitrogen and bacterial diversity (Table 1). Nitrogen is an essential nutrient for microbial growth and plays important roles in controlling microbial diversity and ecosystem productivity (Vitousek et al., 2002; Xia et al., 2008; Sun et al., 2014). In the present study, the nitrogen concentration significantly reduced after the snowfall (Figure 1), which could be due to the snowmelt, which lead nitrogen to be flushed out. Generally, the rate at which chemical compounds are released at snowmelt depends on the way a chemical is incorporated in the snow grain (Pomeroy et al., 2005). Important microbial nutrients, such as NO₃^–, are not incorporated into the snow grain structure and are released easily (Tranter et al., 1986). The positive correlations between nitrogen availability and alpha diversity indices suggest that nitrogen limitation is an important determinant of bacterial diversity. Similar conclusions have been reported in three glaciers on Tibetan Plateau that TN concentrations were significantly correlated with the alpha diversity (Liu et al., 2017).

Consistent with the alpha-diversity indices, the bacterial community changed after snowfall disturbance (Figure 4), again indicating greater environmental filtering effect. The PCoA results showed that the samples on the 26th day were more similar to those before snowfall, while the samples on the 31st day were significantly different. This may be due to a delayed response to disturbance of the microbial community in cryoconite samples. This was also consistent with the increased deterministic processes in bacterial community structure after snowfall (Figure 7). Increased determinism is frequently attributed to the enhanced environmental filtering (Stegen et al., 2012). Our results demonstrated that TN significantly affected the bacterial community (Table 2). This is consistent with the previous findings that nitrogen availability strongly regulates microbial community structure in cryoconite and ice samples from Arctic and Tibetan Plateau glaciers (Christner et al., 2003; Mueller and Pollard, 2004; Edwards et al., 2011; Cameron et al., 2012; Liu et al., 2017). For example, a prior study indicated that nitrite and nitrate concentrations had important effects on the microbial community local scale variations in cryoconite on the Canada Glacier (Antarctic) (Mueller and Pollard, 2004).

Biotic interactions can explain a substantial proportion of the community structure variations (Hacquard et al., 2015; Dang and Lovell, 2016). A previous study of cryoconite from the High Arctic Foxfonna ice cap found that the assembly of the bacterial community of cryoconite was driven principally by biotic filtering (Gokul et al., 2016). Our results indicated that the community network was more complex after the snowfall disturbance, as evidenced by the higher nodes and edge numbers (Table 3). This is likely due to the enhanced environmental filtering, which has been observed in other systems subjected to nitrogen limitation (Wang et al., 2018). Positive correlations were substantial components in the networks before and after the snowfall (Table 3). Positive relationships could be used to infer commensal, mutualistic, obligate syntrophic, and parasitic relationships between bacterial species (Ji et al., 2019). The increased negative correlations after snowfall disturbance may be due to nutrient competition caused by nitrogen limitation (Zhao et al., 2016). It was proposed that the relationship between microorganisms occupying a similar ecological niche could change from commensalism (positive correlations) to competition (negative correlations) when nutrients become limited (Ji et al., 2019). Thus, the nitrogen limitation may influence bacterial diversity and community structure directly and indirectly.

The rare sub-community possessed a higher diversity than the abundant sub-community (Supplementary Figure 5). Thus, the rare taxa are the reservoir of genetic seeds to maintain community diversity and stability (Lennon and Jones, 2011; Rocca et al., 2020), providing biological buffering capacity to withstand environmental changes (Yachi and Loreau, 1999). Our results support prior studies that the rare bacterial taxa have an opportunity to be activated to maintain the stability of the bacterial community under nutrients limitations (Jousset et al., 2017). The snowfall disturbance leads to a decrease in nitrogen content and there is a decrease in temperature at the same time. After the snowfall disturbance, some rare species were adapted to the cold and oligotrophic environment and became abundant species, contributing 22% to the overall community. Four ASVs belonged to Flavobacterium (ASV_4, OTU_26, OTU_145, and OTU_146), one ASV belonged to Psychrobacter (ASV_9) belong to the rare species in samples before snowfall became abundant after snowfall disturbance (Supplementary Table 3). The genus Flavobacterium is known as exopolysaccharide (EPS) producing psychrophiles bacterium. EPS can substantially improve the tolerance of Flavobacterium from freeze–thaw cycles to survive in cold regions. When exposed to a cold environment suddenly, a number of physiological changes occur in the cells, and a set of small molecules are synthesized. The genes encoding these cold-shock-inducible proteins, such as cold shock protein (cspA), have been identified in Flavobacterium and Psychrobacter strains (Liu et al., 2019). They have been frequently reported in oligotrophic Antarctic freshwater systems (Michaud et al., 2012), Antarctic cryoconite holes from Canada Glacier (Christner et al., 2003), and Svalbard soils (Hansen et al., 2007), thus suggesting their bloom in oligotrophic and cold conditions after the snowfall.