Temporal variation of bacterial community and nutrients in Tibetan glacier snowpack

The Tibetan Plateau harbors the largest number of glaciers outside the polar regions, which are the source of several major rivers in Asia. These glaciers are also major sources of nutrients for downstream ecosystems, while there is a little amount of data available on the nutrient transformation processes on the glacier surface. Here, we monitored the carbon and nitrogen concentration changes in a snowpit following a snowfall in the Dunde Glacier of the Tibetan Plateau. The association of carbon and nitrogen changes with bacterial community dynamics was investigated in the surface and subsurface snow (depth at 0–15 and 15–30 cm, respectively) during a 9 d period. Our results revealed rapid temporal changes in nitrogen (including nitrate and ammonium) and bacterial communities in both surface and subsurface snow. Nitrate and ammonium concentrations increased from 0.44 to 1.15 mg L⁻¹ and 0.18 to 0.24 mg L⁻¹ in the surface snow and decreased from 3.81 to 1.04 and 0.53 to 0.25 mg L⁻¹ in the subsurface snow over time. Therefore, we suggest that the surface snow is not nitrogen-limited, while the subsurface snow is associated with nitrogen consumption processes and is nitrogen-limited. The nitrate concentration co-varied with bacterial diversity, community structure, and the predicted nitrogen fixation and nitrogen assimilation/denitrification-related genes (narG), suggesting nitrogen could mediate bacterial community changes. The nitrogen limitation and enriched denitrification-related genes in subsurface snow suggested stronger environmental and biotic filtering than those in surface snow, which may explain the lower bacterial diversity, more pronounced community temporal changes, and stronger biotic interactions. Collectively, these findings advance our understanding of bacterial community variations and bacterial interactions after snow deposition and provide a possible biological explanation for nitrogen dynamics in snow.

1 Introduction

The Tibetan Plateau is the world's third-largest ice reservoir after those in Antarctica and Greenland (Qiu, 2012). These glaciers are the source of several large rivers in Asia, such as the Yellow, Yangtze, Mekong, Salween, Brahmaputra, and Indus rivers (Immerzeel et al., 2010). Glaciers are major sources of nutrients (carbon and nitrogen) for the downstream ecosystems (Singer et al., 2012; Hood et al., 2015; Liu et al., 2021). It has been estimated that 80 Gg of dissolved organic carbon and 27–43 Gg of nitrogen are exported from the Greenland Ice Sheet (Bhatia et al., 2013; Wadham et al., 2016). These nutrients are subjected to complex accumulation and transformation processes in the glacier snow before being released into downstream ecosystems, and microorganisms are the drivers of these processes (Anesio and Laybourn-Parry, 2012; Hell et al., 2013; Hodson et al., 2008). Several studies on snowpacks revealed vital knowledge of the nutrient and microbial community dynamics in the Arctic (Hell et al., 2013; Larose et al., 2013a, 2013b; Maccario et al., 2014, 2019), Antarctic (Antony et al., 2016), and Alps (Lazzaro et al., 2015). However, such knowledge is rarely available for the Tibetan Plateau, constraining our understanding of the nutrient accumulation, transformation, and release processes, which is urgently needed under the enhanced warming and glacier retreat on the Tibetan Plateau.

Autochthonous (microbial origin) and allochthonous (wet and dry atmospheric depositions) are the major sources of nutrients in supraglacial snow, and the contribution of allochthonous sources was much greater in Arctic glaciers (Larose et al., 2013a). Microorganisms are highly involved in the transformation of both autochthonous and allochthonous nutrients. Several studies investigated the dynamics of nutrient and bacterial changes in supraglacial snow during the ablation period. Larose et al. (2013a) revealed that the form of nitrogen varied as a function of time in supraglacial snow during a 2-month field study in Svalbard, and fluctuations in microbial community structure were reported with the relative abundance of fungi and bacteria (such as Bacteroidetes and Proteobacteria) increasing and decreasing, respectively. Seasonal shifts in snowpack bacterial communities were reported in the mountain snow in Japan, where rapid microbial growth was observed with increasing snow temperature and meltwater content (Segawa et al., 2005). However, the results of these studies are likely the consequence of several precipitation events due to the long study period. During precipitation, a new snow layer forms above the previous ones, which is responsible for the stratified snowpack structure. These different snow layers have distinct physical and chemical characteristics, and their age also differs substantially (Lazzaro et al., 2015). Thus, while the microbial process across the aged snowpack can be complex, focusing on supraglacial snow from a single snowfall event could provide unique insights into the bacterial and nutrient dynamics. For instance, Hell et al. (2013) reported bacterial community structure changes during the ablation period across 5 d in the high Arctic, but the bacterial and nutrient dynamics during the snow accumulation period remain elusive.

Surface and subsurface snow typically harbors distinct bacterial community structures (Xiang et al., 2009; Møller et al., 2013; Carey et al., 2016). For example, algae (chloroplasts), Proteobacteria, Bacteroidetes, and Cyanobacteria were more abundant in surface snow, while Firmicutes and Fusobacteria were more abundant in the deeper snow layer (Møller et al., 2013). A previous study had proposed that nitrogen availability could also be a driver of microbial community structure and function in snow (Larose et al., 2013b), in which the NO${}_{\mathrm{3}}^{-}$ and NH${}_{\mathrm{4}}^{+}$ concentrations drove the community composition in Ny-Ålesund snowpack. A dissolved inorganic nitrogen addition experiment also showed a clear community response with the bacterial abundance being elevated and genera richness declining in the final time point compared to the initial time point, suggesting potential specialization of heterotrophic communities (Holland et al., 2020).

Differences in physicochemical conditions can also indirectly influence bacterial community structure through impacts on the types of biotic interactions (Friedman and Gore, 2017; Khan et al., 2018; Bergk Pinto et al., 2019). For example, the addition of organic carbon shifted bacterial interactions from collaboration to competition in Arctic snow (Bergk Pinto et al., 2019), with complex organic carbon degradation and mineralization requiring intensive microbial collaborations (Krug et al., 2020), which are particularly important for oligotrophic environments, such as glaciers. Collaboration is also known to be essential to biological processes such as ammonia oxidation and denitrification, in which various organisms carry out different steps of these processes (Henry et al., 2005; Madsen, 2011; Yuan et al., 2021). These changes in interactions and network complexity can favor or disadvantage certain bacterial groups, thereby changing the bacterial community structure (i.e., biofiltering).

Several studies have investigated the nutrient and bacterial community changes in supraglacial snow across the winter (Brooks et al., 1998; Liu et al., 2006), but the bacterial and nutrient dynamics of freshly fallen snow have been largely overlooked. These short temporal changes will influence the following post-depositional processes after it is buried by the next snowfall and will ultimately determine the physicochemical properties of the stratified snow in the following year. In the present study, we investigated the bacterial community and snow physiochemical property changes in the surface and subsurface supraglacial snow during a 9 d period after a single snowfall event at the Dunde Glacier on the northeastern Tibetan Plateau. We aimed to answer the following key questions: (1) do the bacterial community and nutrients change at a short temporal scale; (2) do the bacterial communities in different snow layers exhibit similar community temporal changes; and (3) are the temporal changes in the surface and subsurface snow related to environmental filtering, biotic interactions, or both?

2 Materials and methods

2.1 Site description and sample collection

Snow samples were collected from the ablation zone at the Dunde Glacier (38^∘06^′ N, 96^∘24^′ E; 5325 m above the sea-level) during October and November 2016 (Fig. S1 in the Supplement). Dunde Glacier is located in the Qilian mountain region on the northeastern Tibetan Plateau, and it is continuously monitored by the Institute of Tibetan Plateau Research, Chinese Academic of Sciences. No supraglacial snow was observed on the glacier surface on 10 October when we first arrived at the camp. Snowfall started on 18 October and ended on 23 October. Sampling was conducted over a 9 d period after the snowfall stopped on a small, flat 5 m × 3 m area to reduce the impact of sample heterogeneity due to spatial variations. Snow samples were collected on 24, 25, 26, 27, and 29 October and 2 November (which are referred to as day 1, 2, 3, 4, 6, and 9) until the next snowfall started. This enabled us to monitor the succession of bacterial communities and the chemical changes in snow through time after deposition. The ambient air temperature during the sampling period averaged −8 ^∘C (data available through the European Centre for Medium-Range Weather Forecasts; Fig. S2), and no snow melting was observed over the 9 d period.

On each day, three snow pits were randomly dug within the 5 m × 3 m area, and any two snow pits were 30–50 cm apart. Each snow pit was approximately 30 cm deep, and the snow was further divided equally into the surface and subsurface layers (approximately 15 cm deep for each layer) to get enough snow for DNA extraction, according to Carey et al. (2016). For each snow pit, the top 1 cm in contact with the air was removed using a sterile spoon to avoid contamination, and then surface and subsurface snow were collected using a sterilized Teflon shovel into 3 L sterile sampling bags separately. Approximately 100 mL were used for physicochemical analyses, whereas the rest was used for DNA extraction. A total of 36 samples were collected. Tyvek bodysuits and latex gloves were worn during the entire sampling process to minimize the potential for contamination, and gloves were worn during all subsequent handling of samples. Samples were kept frozen during the transportation to the laboratory and stored at −20 ^∘C until analysis.

2.2 Environmental characterization of snow

The 100 mL snow sample allocated for physicochemical analysis was melted at room temperature for 3 h before being analyzed. For dissolved organic carbon (DOC) and major ion measurements, 100 mL of snow meltwater was syringe-filtered through a 0.45 µm polytetrafluoroethylene (PTFE) membrane filter (Macherey–Nagel) into 20 mL glass bottles. The membrane was pre-treated with 1 % HCl, deionized water rinsed, and 450 ^∘C combusted for >3 h to remove any potential carbon and nitrogen on the membrane, and the initial 10 mL of the filtrate was discarded before collecting the sample for analysis to eliminate any residual compound on the membrane. The DOC concentrations were measured with a TOC-VCPH analyzer (Shimadzu Corp., Japan). Major ions (NH${}_{\mathrm{4}}^{+}$, NO${}_{\mathrm{3}}^{-}$, Na⁺, K⁺, and SO${}_{\mathrm{4}}^{\mathrm{2}-}$) were analyzed using a Thermo-Fisher ICS-900 (ion chromatography system) as described previously (Rice et al., 2012). The precision and accuracy of the TOC-VCPH analyzer were both <3 %, and the limit of detection was 0.05 mg L⁻¹. The precision and accuracy of the ICS-900 were <5 % and 0.1 mg L⁻¹, and the limit of detection was 0.01 mg L⁻¹ (Fig. S3).

2.3 DNA extraction

For assessing the bacterial community composition, snow samples (3 L) were melted at 4 ^∘C overnight and filtered onto a sterile 0.22 µm polycarbonate membrane (Millipore, USA) with a vacuum pump (Ntengwe, 2005). Bacterial community DNA was extracted from the biomass retained on the filters using a FastDNA^® SPIN Kit for Soil (MP Biomedicals, Santa Ana, CA, USA) according to the manufacturer's instructions. DNA extraction with no sample added was performed in parallel and used as a negative control.

The raw DNA was checked by electrophoresis in 1 % ($w/v$) agarose gel and purified from the gel using an Agarose Gel DNA purification kit (Takara, Japan). The concentration and purity of the DNA extracts were measured using a NanoDrop 1000 spectrophotometer (Thermo Scientific, Wilmington, DE, USA). The extracted DNA was stored at −80 ^∘C until amplification.

2.4 Bacterial 16S rRNA amplification and Illumina MiSeq sequencing

In total, 36 DNA samples and one negative control were subjected to amplicon sequencing. Universal primers 515F (5^′-GTGCCAGCMGCCGCGGTAA-3^′) and 806R (5^′-GGACTACHVGGGTWTCTAAT-3^′) (Caporaso et al., 2012), with 12 nt unique barcodes, were used to amplify the V4 hyper-variable regions of the bacterial 16S rRNA gene. Polymerase chain reaction (PCR) was performed under the following conditions: 94 ^∘C for 5 min, 30 cycles of 94 ^∘C for 30 s, 52 ^∘C for 30 s, and 72 ^∘C for 30 s, followed by a final cycle of 10 min at 72 ^∘C. Each PCR reaction contained 12.5 µL 2× Premix Taq DNA polymerase (Takara Biotechnology and Dalian Co. Ltd., China), 1 µL primer (0.4 µM final concentration), and 8.5 µL nuclease-free water, as well as 2 µL DNA template (20 ng µL⁻¹) or 2 µL sterile water for the PCR negative controls. PCR products were confirmed using agarose gel electrophoresis, and no PCR band was detected in the PCR negative controls. To minimize PCR batch-to-batch variations and maximize the quantity of PCR product, triplicate PCR reactions were performed for each sample, and PCR products were pooled for purification using the OMEGA Gel Extraction Kit (Omega Bio-Tek, Norcross, GA, USA) following electrophoresis. PCR products from different samples were pooled in equal molar amounts and then used for 2×250 bp paired-end sequencing on a MiSeq machine (Illumina, San Diego, CA).

2.5 Processing of Illumina sequencing data

MiSeq sequence data were processed using the QIIME 2 pipeline version 2018.8 (Bolyen et al., 2019), following the recommended procedures (https://docs.qiime2.org/2022.2/, last access: 4 March 2020) and using the plugin demux to visualize interactive quality diagrams and check read quality. Plugin DADA2 (Callahan et al., 2016) was applied to remove primers, truncate poor-quality bases, conduct de-replication, identify chimeras, and merge paired-end reads. Commands included in the feature table (McDonald et al., 2012) generated the summary statistics of sequences related to the samples. Further, we trained a naïve Bayes classifier with the feature-classifier plugin using the 16S rRNA gene database at 99 % similarity of the SILVA 132 QIIME release and based on the 515F–806R primer pair as used for the PCR. Finally, the taxa plugin was used to filter mitochondrial and chloroplast sequences, as well as to generate absolute read count tables of all taxa for each sample. Data were analyzed at the level of amplicon sequence variant (ASV), in which ASVs are delineated by 100 % sequence identity (Callahan et al., 2017).

After removing singletons, a total of 1 685 186 high-quality reads were obtained, representing 9178 ASVs. Before statistical analysis, the dataset was rarefied to 45 000 reads per sample, which is the lowest read count among samples. Rarefaction curves reached an asymptote before the subsampling, which confirmed that this depth was sufficient to detect the diversity present (Fig. S4).

2.6 Network analysis

The ASV–ASV associations within the surface and subsurface bacterial communities were explored using Molecular Ecological Network Analysis Pipeline (http://ieg4.rccc.ou.edu/mena, last access: 10 May 2021) (Deng et al., 2012). The ASVs that occurred in at least 50 % of the samples from the surface or subsurface group were selected to construct the network. Spearman's rank correlation coefficient (ρ) was calculated to reflect the strength of association between species. The false discovery rates (q values) were calculated from the observed p-value distribution. The resulting correlation matrix was analyzed with the random matrix theory (RMT)-based network approach to determine the correlation threshold for network construction, and the same threshold was used for both the surface and subsurface network, so the topological properties of the surface and subsurface networks are comparable.

2.7 Statistical analysis

Shannon–Wiener and Chao1 indices, which were used to estimate the species richness in the snow community, were calculated using the diversity function in the R package vegan (Oksanen et al., 2010). Functional profiling of bacterial taxa was carried out using the package Tax4Fun2 in R (Wemheuer et al., 2020). While the application of functional profiles predicted from 16S rRNA gene-based community composition data is limited by the functional information available in databases, we present these data as one possible interpretation of the patterns detected and note that the Tax4Fun2 package performed well compared to older widely used programs (Wemheuer et al., 2020). The pairwise Wilcoxon rank-sum test was used to compare the depth-horizon differences in environmental variables, alpha diversity, and the relative abundance of taxonomic groups at the phylum level. Linear regression modeling was implemented in R using the lm function to estimate the trend of environmental characteristics, alpha diversity, and microbial community composition changes. Multiple linear regression analysis was performed to determine the contribution and significance of the environmental characteristics to the alpha diversity using the lm function in R. We use the stepwise Akaike information criterion (AIC) method for variable selection by the step function in R. The best model was chosen based on the lowest AIC value (Wagenmakers and Farrell, 2004). The bacterial community structure was subjected to principal coordinate analysis (PCoA) carried out using the pcoa function of the ape package in R. The significance of dissimilarity of community composition among samples 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 statistically significant. The Mantel test based on Spearman's rank correlations was performed using the bacterial dissimilarity and environmental dissimilarity matrices, calculated based on the Bray–Curtis distance metrics and Euclidean distance metrics in the vegan R package, respectively. 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. All environmental variables were normalized before the calculation. All statistical analyses were executed in R version 3.4.3 (R Core Team, 2017).

Figure 1 The pattern of environmental factor changes in the surface and subsurface snow layers. (a) Environmental factor comparisons in the surface and subsurface snow layers. Each dot represents an individual sample. Significantly higher concentrations of NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$, K⁺, and SO${}_{\mathrm{4}}^{\mathrm{2}-}$ were observed in the subsurface layer based on Wilcoxon rank-sum test. (b) Temporal changes in environmental factors in the surface and subsurface layers. The solid and dashed lines indicate significant and non-significant temporal changes, respectively. The concentration of NO${}_{\mathrm{3}}^{-}$ and NH${}_{\mathrm{4}}^{+}$ in the surface layer significantly increased with time, while the concentration of NO${}_{\mathrm{3}}^{-}$, and NH${}_{\mathrm{4}}^{+}$ in the subsurface layer, significantly decreased with time. Significance is based on linear regression. Grey shading indicates the 95 % confidence interval of regression.

3 Results

3.1 Environmental characteristics of the snowpack

The concentrations of NO${}_{\mathrm{3}}^{-}$ and NH${}_{\mathrm{4}}^{+}$ ranged from 0.44 to 5.09 and 0.17 to 0.62 mg L⁻¹, respectively (Fig. 1a, Table S1 in the Supplement), and they were both significantly higher in the subsurface than in the surface snow (Wilcoxon rank-sum test: all p<0.001; Fig. 1a). K⁺ and SO${}_{\mathrm{4}}^{\mathrm{2}-}$ ions in the subsurface snow were also significantly higher (0.29±0.13 and 6.09±3.18 mg L⁻¹, respectively) than those in the surface snow (0.12±0.08 and 3.71±1.64 mg L⁻¹; Wilcoxon rank-sum test: p<0.001 and p=0.015, respectively). The concentrations of DOC ranged from 0.46 to 5.89 mg L⁻¹ and exhibited no significant difference between the surface and subsurface snow (Wilcoxon rank-sum test: p=0.310). The concentrations of Na⁺ ions ranged from 0.35 to 7.34 mg L⁻¹, and there was no significant difference between the surface and subsurface snow (Wilcoxon rank-sum test: p=0.079). The concentration of NO${}_{\mathrm{3}}^{-}$ and NH${}_{\mathrm{4}}^{+}$ ions in the surface snow exhibited a weak but significantly positive association with time ($F_{\mathrm{1},\mathrm{16}}=\mathrm{5.97}$, p=0.027, and R²=0.27 and $F_{\mathrm{1},\mathrm{16}}=\mathrm{8.58}$, p=0.010, and R²=0.35, respectively; Fig. 1b). On the other hand, stronger negative associations were found between inorganic nitrogen and time in the subsurface snow ($F_{\mathrm{1},\mathrm{16}}=\mathrm{40.66}$, p<0.001, and R²=0.72 and $F_{\mathrm{1},\mathrm{16}}=\mathrm{50.74}$, p<0.001, and R²=0.76, respectively). Other environmental parameters exhibited no significant changes with time.

Figure 2 Taxonomic composition of bacterial community in snow. Only dominant phyla are presented (relative abundance >1 %). The snow community are dominated by Alphaproteobacteria, Actinobacteria, Cyanobacteria, Gammaproteobacteria, Bacteroidetes, Firmicutes, Chloroflexi, Gemmatimonadetes, Planctomycetes, Acidobacteria, Deltaproteobacteria, and Deinococcus-Thermus.

3.2 Diversity and composition of bacterial community from the snowpack

The surface and subsurface snow were both dominated by Alphaproteobacteria, Actinobacteria, Cyanobacteria, Gammaproteobacteria, Bacteroidetes, Firmicutes, Chloroflexi, Gemmatimonadetes, Planctomycetes, Acidobacteria, Deltaproteobacteria, and Deinococcus-Thermus (Fig. 2). The relative abundance of most of these phyla was not significantly different in the two snow layers except the Gemmatimonadetes, Planctomycetes, and Acidobacteria, which were significantly more abundant in the surface layer than in the subsurface layer (all p<0.05, Wilcoxon rank-sum test; Fig. S5). In the surface layer, weak but significant negative trends were observed between the relative abundances and ASV number of Alphaproteobacteria, Gammaproteobacteria, and Firmicutes and time ($F_{\mathrm{1},\mathrm{16}}=\mathrm{6.97}$, p=0.018, and R²=0.30; $F_{\mathrm{1},\mathrm{16}}=\mathrm{23.8}$, p<0.001, and R²=0.60; and $F_{\mathrm{1},\mathrm{16}}=\mathrm{22.28}$, p<0.001, and R²=0.58 in relative abundance; $F_{\mathrm{1},\mathrm{16}}=\mathrm{7.56}$, p=0.014, and R²=0.32; $F_{\mathrm{1},\mathrm{16}}=\mathrm{27.12}$, p<0.001, and R²=0.63; and $F_{\mathrm{1},\mathrm{16}}=\mathrm{16.68}$, p=0.001, and R²=0.51 in ASV number, respectively), while weak positive correlations were observed between the relative abundances and ASV number of Cyanobacteria and Deinococcus-Thermus and time ($F_{\mathrm{1},\mathrm{16}}=\mathrm{6.94}$, p=0.018, and R²=0.30 and $F_{\mathrm{1},\mathrm{16}}=\mathrm{13.10}$, p=0.002, and R²=0.45 in relative abundance; $F_{\mathrm{1},\mathrm{16}}=\mathrm{3.42}$, p=0.083, and R²=0.18 and $F_{\mathrm{1},\mathrm{16}}=\mathrm{4.07}$, p=0.061, and R²=0.20 in ASV number, respectively; Figs. S6 and S7 in the Supplement). Relative to the surface snow, the subsurface layer had a stronger negative correlation between the relative abundance and ASV number of Alphaproteobacteria and Firmicutes and time ($F_{\mathrm{1},\mathrm{16}}=\mathrm{15.17}$, p=0.001, and R²=0.49 and $F_{\mathrm{1},\mathrm{16}}=\mathrm{15.43}$, p=0.001, and R²=0.49 in relative abundance; $F_{\mathrm{1},\mathrm{16}}=\mathrm{18.98}$, p=0.083, and R²=0.54 and $F_{\mathrm{1},\mathrm{16}}=\mathrm{15.17}$, p=0.001, and R²=0.53 in ASV number, respectively; Figs. S6 and S7), while weak correlations were observed between the relative abundance and ASV number of Cyanobacteria and Chloroflexi and time ($F_{\mathrm{1},\mathrm{16}}=\mathrm{5.62}$, p=0.031, and R²=0.26 and $F_{\mathrm{1},\mathrm{16}}=\mathrm{12.81}$, p=0.003, and R²=0.44 in relative abundance; $F_{\mathrm{1},\mathrm{16}}=\mathrm{5.34}$, p=0.034, and R²=0.25 and $F_{\mathrm{1},\mathrm{16}}=\mathrm{14.49}$, p=0.002, and R²=0.47 in ASV number, respectively).

Figure 3 Bacterial alpha diversity in snow layers. (a) Bacterial alpha diversity comparison between the surface and subsurface layers. Each dot represents an individual sample. For both Shannon and Chao1 indices, no significant difference was observed between the surface and subsurface snow layers. The comparison is based on Wilcoxon rank-sum test. (b) Temporal changes in the alpha diversity indices in the surface and subsurface snow layers. For the surface layer, no significant correlation was observed, while both Shannon and Chao1 showed a significantly reduction with time in the subsurface layer. Significance is based on linear regression. Grey shading indicates the 95 % confidence interval of regression.

Figure 4 The influence of environmental factors on bacterial diversity. Correlations of Shannon (a, c) and Chao1 (b, d) diversity indices with environmental factors in the surface and subsurface layers. Each dot represents an individual sample. The solid and dashed lines indicate significant and nonsignificant changes, respectively. Significance is based on linear regression. Grey shading indicates the 95 % confidence interval of regression.

The bacterial Shannon and Chao1 indices in the surface snow were 5.61±0.39 and 744±199, respectively, and were not significantly different from those in the subsurface layer (5.52±0.68 and 705±269, respectively) (p=0.81 and 0.57, respectively) (Fig. 3a). In the surface snow, the Shannon and Chao1 indices were similar across the 9 d ($F_{\mathrm{1},\mathrm{16}}=\mathrm{0.37}$, p=0.553, and R²=0.02 and $F_{\mathrm{1},\mathrm{16}}=\mathrm{0.01}$, p=0.939, and R²=0.001, respectively; Fig. 3b). Besides, weak positive associations of Shannon and Chao1 indices with the DOC and sodium ions were detected ($F_{\mathrm{1},\mathrm{16}}=\mathrm{4.90}$, p=0.042, and R²=0.23 and $F_{\mathrm{1},\mathrm{16}}=\mathrm{4.91}$, p=0.042, and R²=0.24, respectively; Fig. 4a and b). In contrast, although weak, significant negative correlations were observed in both Shannon and Chao1 indices with time in the subsurface snow ($F_{\mathrm{1},\mathrm{16}}=\mathrm{12.33}$, p=0.003, and R²=0.44 and $F_{\mathrm{1},\mathrm{16}}=\mathrm{8.73}$, p=0.009, and R²=0.35, respectively). Weak but significant positive associations of Shannon and Chao1 indices with the concentrations of NO${}_{\mathrm{3}}^{-}$ and NH${}_{\mathrm{4}}^{+}$ were detected (Shannon diversity: $F_{\mathrm{1},\mathrm{16}}=\mathrm{9.13}$, p=0.008, and R²=0.36 and $F_{\mathrm{1},\mathrm{16}}=\mathrm{5.17}$, p=0.037, and R²=0.24, respectively; Chao1 index: $F_{\mathrm{1},\mathrm{16}}=\mathrm{8.60}$, p=0.009, and R²=0.36 and $F_{\mathrm{1},\mathrm{16}}=\mathrm{5.32}$, p=0.035, and R²=0.25, respectively; Fig. 4c and d). This is consistent with the multiple linear regression results, which consistently identified the concentrations of NO${}_{\mathrm{3}}^{-}$ and NH${}_{\mathrm{4}}^{+}$ as the significant determinants of bacterial Shannon diversity in the subsurface layer (Table S2).

Figure 5 Principal coordinate analysis (PCoA) of microbial communities in the surface and subsurface snow. (a) Bray–Curtis distance-based PCoA ordination plot. The microbial community structures of the surface and subsurface snows are significantly different (PERMANOVA, p<0.001). (b) Pairwise regression analysis between PCoA scores and sampling time. The solid and dashed lines indicate significant and insignificant changes (based on linear regression), respectively. The PCoA1 scores for the bacterial community in the surface layer exhibit no significant correlation with time, while the PCoA2 scores significantly correlated with time. The PCoA1 and PCoA2 are both significantly correlated with time in the subsurface layer. Grey shading indicates the 95 % confidence interval of regression.

3.3 Bacterial community structure and functional genes

The bacterial community structure at the ASV level significantly differed in the surface and subsurface snow (PERMANOVA, F=2.78, p<0.001; Fig. 5a), as well as among the different sampling times (PERMANOVA, F=3.31 and p<0.001 and F=2.17 and p<0.001, respectively). Additionally, a significant interactive effect was detected between the depth and time (PERMANOVA, F=2.68, p<0.001), indicating that the depth influenced the temporal pattern of bacterial community structure changes. Specifically, only the second principal coordinate (PCoA2) values of the surface snow significantly varied with time ($F_{\mathrm{1},\mathrm{16}}=\mathrm{141.8}$, p<0.001, R²=0.89; Fig. 5b), while the PCoA1 values of the surface snow did not ($F_{\mathrm{1},\mathrm{16}}=\mathrm{0.04}$, p=0.840, R²=0.003; Fig. 5b). Furthermore, PCoA1 and PCoA2 of the surface snow exhibited no significant correlation with the measured environmental factors (all p>0.05; Figs. S8 and S9). In comparison, both PCoA1 and PCoA2 values of the subsurface, albeit weakly, co-varied with time ($F_{\mathrm{1},\mathrm{16}}=\mathrm{6.35}$, p=0.023, and R²=0.28 and $F_{\mathrm{1},\mathrm{16}}=\mathrm{8.38}$, p=0.011, and R²=0.34, respectively; Fig. 5b), while the PCoA2 also demonstrated a significant association with nitrate, ammonium, potassium, sulfate, and DOC concentrations (all p<0.05; Fig. S9).

Normalized stochasticity ratio (NST) was used to examine the relative contributions of stochasticity and determinism in shaping bacterial communities. The average NST values were 74 % and 46 % in the surface and subsurface snow layers, and the contribution of stochasticity was significantly higher in the surface than in the subsurface layers (p<0.001; Fig. S10).

Table 1 Results of Mantel test showing the relationships between bacterial community composition and environmental factors in the surface and subsurface snow. Significant correlations are in bold.

Mantel tests were performed to evaluate the effects of environmental factors on bacterial community structure for each layer. No significant correlation was identified between the measured environmental factors and the bacterial community structure in the surface snow. However, weak positive associations were apparent in the subsurface snow with the concentrations of NO${}_{\mathrm{3}}^{-}$ and NH${}_{\mathrm{4}}^{+}$ (p=0.005 and 0.01, respectively) (Table 1). The relative abundance of nitrogen-cycling-associated functional genes was predicted in the surface and subsurface snow. The relative abundance of the nitrogen-fixation marker gene (nifH) positively associated with time in the surface layer, while no clear pattern was observed in the subsurface layer ($F_{\mathrm{1},\mathrm{16}}=\mathrm{7.76}$, p=0.013, and R²=0.33 and $F_{\mathrm{1},\mathrm{16}}=\mathrm{0.57}$, p=0.461, and R²=0.01, respectively; Fig. S11). The relative abundance of the narG gene, which is involved in the nitrate reduction and denitrification process, exhibited negative and positive associations with time in the surface and subsurface, respectively ($F_{\mathrm{1},\mathrm{16}}=\mathrm{4.69}$, p=0.046, and R²=0.23 and $F_{\mathrm{1},\mathrm{16}}=\mathrm{11.24}$, p=0.004, and R²=0.41, respectively). The nirK gene, which is also involved in the denitrification process, decreased with time in the surface layer, while no significant change was detected in the subsurface layer ($F_{\mathrm{1},\mathrm{16}}=\mathrm{10.39}$, p=0.005, and R²=0.39 and $F_{\mathrm{1},\mathrm{16}}=\mathrm{1.98}$, p=0.179, and R²=0.05, respectively).

Figure 6 Bacterial co-occurrence networks for the surface and subsurface layer communities. Each node represents a bacterial amplicon sequence variant (ASV). The solid red lines represent positive correlations, and the solid blue lines represent negative correlations. Nodes are colored by taxonomy at the phylum level. The subsurface community networks are more complex with a higher positive-to-total correlation ratio.

Table 2 Topological properties of the empirical networks for the surface and subsurface bacterial communities.

3.4 Interspecies interactions at the surface and subsurface layers

Co-occurrence networks were constructed for the surface and subsurface bacterial communities to infer the biotic interactions among species (Fig. 6). The surface network comprised a higher number of nodes (each indicating one ASV, node number =197) but a lower number of edges (each indicating a significant association between two ASVs; edge number =436) than the subsurface network (node number =140 and edge number =523, Table 2). The network in the subsurface snow, relative to surface snow, demonstrated a higher number of edges per node (3.73 and 2.21, respectively), higher average connectivity (avgK; 7.57 and 4.43, respectively), and lower average path distance (GD; 4.72 and 5.51, respectively), which indicate a substantially more complex network topology. Both networks were dominated by positive (co-presence) relationships, and the subsurface network exhibited a higher positive-to-total interaction ratio (95 %) than the surface network (83 %).

Modularity, average clustering coefficient (avgCC), and graph density of the surface and subsurface bacterial community networks were all higher than those of random networks (Table S3), indicating that snowpack bacterial networks showed non-random assemblage and exhibited modular structures. The subsurface networks showed higher values of avgCC (0.39), transitivity (0.49), and connectedness (0.86) than the surface bacterial community network (0.31, 0.45, and 0.71, respectively), indicating a greater degree of connectivity (Table 2).

4 Discussion

4.1 Rapid shifts of bacterial community structure across a short temporal scale

The surface and subsurface snow was dominated by Alphaproteobacteria, Actinobacteria, Cyanobacteria, Gammaproteobacteria, and Bacteroidetes (Fig. 2). Despite differences in sampling season, the bacterial taxa detected were consistent with previous studies on snow in the Arctic and Antarctic (Larose et al., 2010; Carpenter et al., 2000; Amato et al., 2007; Lopatina et al., 2013; Møller et al., 2013). Bacterial richness and diversity exhibited little change throughout the 9 d in the surface snow layer, while they exhibited a reduction trend in the subsurface snow layer (Fig. 3b). This indicates that the microbiome in the subsurface snow may be subjected to greater environmental filtering than those in the surface snow (Xiang et al., 2009). Among all environmental factors measured, nitrate and ammonium were the only measured environmental factors that changed across the 9 d. The nitrate and ammonium concentrations in the subsurface snow both exhibited an R² value of greater than 0.7 and reduced with time, therefore indicating a consumption process (Fig. 1b). Despite the R² value being weak, both nitrate and ammonium concentrations co-varied with bacteria richness and diversity in subsurface snow, which was not observed in the surface snow (Fig. 4). Furthermore, multiple linear regression analyses also identified nitrate and ammonium to be the dominant driver of bacterial Shannon diversity in the subsurface snow (Table S2). Thus, these results suggest that nitrate and ammonium could play a more important role in influencing bacterial diversity in subsurface snow than that in surface snow. Nitrogen is an essential nutrient for microbial growth, and it plays an important role in controlling microbial diversity and ecosystem productivity (Vitousek et al., 2002; Xia et al., 2008; Sun et al., 2014). The positive associations between nitrogen concentration and alpha diversity indices have been typically inferred as nitrogen limitation (Telling et al., 2011). Thus, these results hint that nitrogen limitation could occur in subsurface snow and influence bacteria diversity. In comparison, the surface layer is unlikely to be subjected to nitrogen limitation, and the nitrogen in the surface snow slightly increased. This is consistent with previous studies on the Greenland ice sheet, where nitrate additions to surface ice did not alter the cryoconite community cell abundance and 16S rRNA gene-based community composition (Cameron et al., 2017).

The bacterial community structure also exhibited temporal changes in the subsurface layer. Furthermore, associations between nitrogen and the microbial community structure were observed to a certain degree (Table 1 and Fig. 5), again indicating some level of environmental filtering (Kim et al., 2016). This is consistent with the finding in the Arctic that nitrogen influences snow bacterial community composition via regulating algae metabolism (Lutz et al., 2017). This is also consistent with the higher contribution of deterministic processes in the subsurface layer than in the surface layer (Fig. S10). Deterministic processes could be due to environmental filtering or biotic interactions, whereas stochastic processes include dispersal limitation, community drift, and speciation (Stegen et al., 2012). The surface layer could receive nitrogen input through aeolian deposition processes (Björkman et al., 2014), whereas the subsurface snow could only receive limited external microbial and nutrient input through supraglacial meltwater. The latter could be particularly limited during the glacier deposition period when the glacier surface temperature is below 0 ^∘C (Fig. S2).

Our results suggest that both bacteria and snow physiochemical properties experience changes across the 9 d during the snow deposition period for the Tibetan glacier investigated here, and those changes were stronger in the subsurface layer than in the surface layer. Traditionally, supraglacial snow is recognized as a cold oligotrophic environment with a very slow metabolism rate (Quesada and Vincent, 2012; Marshall and Chalmers, 1997), but increasing evidence has suggested that bacterial community changes can occur on a short temporal scale. For example, Hell et al. (2013) reported changes in the dominant bacterial phylum Proteobacteria across 5 d, and active bacterial metabolism has been observed in the Greenland Ice Sheet supraglacial ice (Nicholes et al., 2019). In addition, active bacteria affiliated with Proteobacteria have been identified in the Antarctic (Lopatina et al., 2013) and Arctic (Holland et al., 2020) snow at temperatures below 0 ^∘C, therefore supporting the present study that bacterial community changes in 9 d could be possible. This indicates that supraglacial snow can harbor an active bacterial community, which in turn can have an impact on nutrient transformation.

4.2 Distinct nitrogen-transformation processes in surface and subsurface snow

Both ammonium and nitrate concentrations showed a weak increasing trend with time in the surface snow (Fig. 1). The weak increase in ammonium could be explained by biogenic emissions due to local plant and animal sources (Filippa et al., 2010), while the increase in nitrate has been largely attributed to atmospheric deposition (Björkman et al., 2014). Nitrogen deposition occurs at a rate of 282 kg N km⁻² yr⁻¹ in the region of our investigation (Lü and Tian, 2007), which equals 0.19 mg N for the 0.5 m × 0.5 m area sampled each day (assuming nitrogen deposition occurred evenly across the year). Another potential source of nitrogen input could be the nitrogen fixation process (Telling et al., 2011). Bacteria are the only microorganisms that are capable of fixing atmospheric nitrogen (Bernhard, 2010). Potential nitrogen input from microbial processes is supported by the increase in the nitrogen-fixing Cyanobacteria (Fig. S6) and nifH gene (Fig. S11). Cyanobacteria are known as free-living phototrophs capable of nitrogen fixation, especially in extreme environments (Chrismas et al., 2018; Makhalanyane et al., 2015; Levy-Booth et al., 2014). For example, Cyanobacteria were found to be the main group of potential nitrogen fixers determined by quantitative PCR with three sets of specific nifH primers on the surface of the Greenland Ice Sheet (Telling et al., 2012). The nitrogen fixation rate was not quantified in the present study, but the present study suggests that microbial nitrogen fixation could be an overlooked source of nitrogen in Tibetan glacier snow. Further transcriptomic and nitrogen-isotope analyses may provide additional evidence on the microbial activity in nitrogen fixation.

In contrast with the surface layer, nitrogen concentrations (nitrate and ammonium) significantly decreased in the subsurface snow with time (Fig. 1). A possible explanation for this might be the microbial utilization and photochemical degradation of nitrogen compounds (Björkman et al., 2014). The microbial processes, i.e. nitrate reduction and denitrification process, are evidenced by the increase in the narG gene (Fig. S11) (Telling et al., 2011; Zhang et al., 2020). Alternatively, microorganisms may carry out assimilatory nitrate reduction, which is used to incorporate nitrogen into biomolecules (Larose et al., 2013a; Richardson and Watmough, 1999). The assimilatory process is performed by a range of microorganisms including bacteria, algae, yeasts, and fungi (Huth and Liebs, 1988). Thus, further studies on eukaryotes, including algae, may provide a full understanding of the nitrogen consumption mechanisms in subsurface snow. The denitrification process converts nitrate to N₂ and generates nitrite, nitric oxide (NO), and nitrous oxide (N₂O) intermediates (Kuypers et al., 2018). A previous study detected microbial-specific phylogenetic probes that targeted genera whose members are able to carry out denitrification reactions such as Roseomonas in a snowpack of Spitsbergen Island of Svalbard, Norway (Larose et al., 2013a). Amoroso et al. (2010) also proposed that denitrification can explain the microbial isotopic signature observed in winter snow at Ny-Ålesund. Although the oxygen level in the subsurface snow was not measured, the occurrence of anaerobic denitrification reactions in subsurface snow has been reported in Arctic snowpacks (Larose et al., 2013a). Lastly, photochemical degradation of nitrogen compounds is the most well-known nitrogen degradation pathway, and the release of both NO and NO_x by NO${}_{\mathrm{3}}^{-}$ photolysis on natural snow has been reported in European high Arctic snowpack (Amoroso et al., 2010; Beine et al., 2003). In a snow reactive nitrogen oxide (NO_y) survey in Greenland, NO_y flux was reported to exit snow in 52 out of 112 measurements (Dibb et al., 1998). Further metatranscriptomic analyses targeting the genes associated with nitrogen cycling are required to confirm the distinct nitrogen transformation processes between the surface and subsurface layers.

4.3 Subsurface snow exhibits greater complexity in biotic interactions

Biotic interactions can explain a substantial proportion of the community structure variations (Hacquard et al., 2015; Dang and Lovell, 2016). Our results indicated that the subsurface community network was more complex as evidenced by the higher average connectivity and a shorter path length (GD) compared to the surface community network (Table 2). This is likely due to the enhanced environmental filtering, as has been observed in other systems subjected to environmental stresses (Ji et al., 2019; Wang et al., 2018). A higher ratio of positive-to-total interactions, but lower modularity, was identified in the subsurface snow network (Table 2). In general, higher positive interactions indicate increased microbial cooperation (Ju et al., 2014; Scheffer et al., 2012), whereas a reduction in modularity indicates microbial niche homogenization (Ji et al., 2019). The enhanced biotic associations and cooperation in the subsurface layer may be attributed to the occurrence of denitrification processes, as denitrification is a multi-step process that involves multiple bacterial cohorts to complete the process (Henry et al., 2005; Madsen, 2011; Yuan et al., 2021). The enhanced collaboration and deterministic succession were previously reported in the bacterial community associated with the anoxic decomposition of microcystis biomass (Wu et al., 2020), while cross-feeding was shown to enhance positive interactions among the different members of the community (Borchert et al., 2021).

The path lengths of the subsurface network were lower than that of the surface layer (Table 2). The shorter path length has been proposed to be associated with a higher transfer efficiency of information and materials across the microorganisms in the network (Du et al., 2020) which are required for complex biological processes that require extensive bacterial collaboration, such as denitrification (Yuan et al., 2021). Thus, the short path length is consistent with the dominance of denitrification processes in the subsurface layer. Previous studies have proposed microbial interactions as biotic drivers that impact microbial diversity (Calcagno et al., 2017; Hunt and Ward, 2015). Thus, those microorganisms which are not adapted to the subsurface environment would be excluded from the environment, which provides an alternative explanation for the reduction in diversity (Scheffer et al., 2012; Ziegler et al., 2018; Bergk Pinto et al., 2019).

5 Conclusion

Our results showed the dynamics of nitrogen and the bacterial community in supraglacial snow over 9 d. Inorganic nitrogen was unchanged or slightly increased in the surface snow, while it decreased in subsurface snow. Due to atmospheric nitrogen deposition and potential bacterial nitrogen fixation activities, nitrogen limitation is unlikely to occur in the surface snow. In contrast, nitrogen consumption was inferred in the subsurface snow. Nitrogen is traditionally recognized to be released from the supraglacial environment due to photolysis, whereas this study hints that nitrogen assimilation and denitrification could be alternative routes. Therefore, the increased nitrogen deposition due to anthropogenic activities may enhance the nitrogen consumption in the subsurface snow, which reduces the impact of increased nitrogen discharge on downstream glacier-fed rivers. In summary, our results provide a new perspective on the nutrients and bacterial community dynamics in supraglacial snow of the Tibetan Plateau. Further studies based on metagenome and metatranscriptome can enhance the understanding of bacterial functions.