Introduction

Glaciers are solid reservoirs that store 75% of the freshwater on the Earth’s surface and have an important influence on the global water cycle1. Global warming has increased the rate and extent of the melting of alpine glaciers over the last few decades2. Glaciers in the Tibetan Plateau are rapidly shrinking, remarkably influencing downstream ecosystems and approximately 2 billion people3. Glacial ecosystems embody a series of unique habitats that harbor distinctive organisms and biogeochemical activity4,5. These ecosystems can be categorized into supraglacial, englacial, subglacial, and proglacial ecosystems, depending on their spatial position relative to the glacier6,7,8. These glacial ecosystems are separately located on the surface, inner, bottom, and front of glaciers, but are connected by glacier meltwater9. Microbes in these glacial ecosystems exhibit varying community structures and disperse among ecosystems10,11. The connection of microbes across different ecosystems is found to increase under changing hydrological conditions due to glacier melting7. Therefore, comparing microbial community diversity and composition among different glacial ecosystems will provide a unique perspective for understanding the connections and shifts of glacial ecosystems as well as their impact on adjacent ecosystems under global warming.

Different glacial ecosystems have distinct temperatures, levels of solar radiation, and nutrient content12,13, which shape diverse microbial communities. The supraglacial ecosystem is oxygenated and illuminated by sunlight, predominantly populated by photoautotrophic microorganisms like cyanobacteria and microalgae14, which are actively involved in microbial carbon15,16 and nitrogen cycling17,18. The subglacial ecosystems are characterized by a paucity of light and restricted oxygen availability at the bottom of glaciers19, harboring taxonomically and functionally diverse prokaryotic microorganisms, including chemoautotrophs and heterotrophs20,21. Studies of subglacial ecosystems usually focus on basal ice and sediments that are collected beneath or near glacier termini due to the difficulty in collecting samples under deep ice cover22,23,24. Supraglacial and subglacial ecosystems are potential microbial sources for proglacial ecosystems through glacier meltwater25,26,27,28. Irvine-Fynn et al.29 estimated that 1.8 × 1018 and 8.3 × 1018 microbial cells were delivered to the subglacial ecosystem by supraglacial stream during ablation season in the Greenland Ice Sheet29. Kohler et al.30 demonstrated that subglacial ecosystems exported regionally specific communities and probably reflected differences in dominant hydrological flow paths in Arctic glaciers30. Proglacial ecosystems, including proglacial lakes, runoff, and foreland, are sensitive to glacier melting and receive upstream microbes and nutrients7. Proglacial lakes receive large plumes of sediment-laden meltwater and are dominated by heterotrophic and eventually mixotrophic microbes31. Moreover, proglacial lakes are increasingly expanding due to rapid glacier melt, posing the risk of proglacial lake outburst floods32.

Microbial communities of glacial ecosystems differ between distinct habitats but also connect and interplay by glacial meltwater30,31,33,34,35,36,37. For instance, snow was a potential microbial seed bank of the glacial-fed lake in the Tibetan Plateau36. Subglacial environments have the potential to provide continuous sources of chemical energy (i.e., H2 and oxidized minerals) to support chemolithotrophic microbial primary production in proglacial environments23. However, little work dealt with the difference between supraglacial and subglacial ecosystem’s impact on proglacial ecosystems. Distinguishing the difference in the contributions of supraglacial and subglacial ecosystems to proglacial ecosystems is meaningful to provide a baseline for forecasting the responses of glacial ecosystems under global warming as well as the impact of glacier retreat on downstream ecosystems. In addition, microbial community diversity and community structure are key indicators of community stability38,39. Network properties (e.g., cohesion and robustness) have been used to successfully predict the stability of microbial community40,41,42,43. However, the stability has rarely been quantified in previous glacial microbial community studies, and it is important to understand the stability of different glacial ecosystems under global warming.

In this study, we collected 42 samples from seven habitats of three glacial ecosystems (i.e., supra-, sub-, and proglacial ecosystems) in the Kuoqionggangri Glacier on the Tibetan Plateau, which is experiencing rapid shrinkage. As the supraglacial and subglacial environments are distinct and harbor different microbial communities, while microbes from different sources could converge downstream due to glacier melting26,27,30, we hypothesized that supraglacial and subglacial species would be exported to the proglacial ecosystem and exhibit different contributions. Testing this hypothesis will help to improve our understanding of the impact of glacier melting on downstream ecosystems. We compared the diversity, taxonomic composition, structure, and stability of microbial communities in different glacial habitats and ecosystems by 16S rRNA gene sequencing. Here, we aim to elucidate the similarities and differences in microbial communities between different glacial ecosystems and to distinguish the contributions of supraglacial and subglacial ecosystems to the proglacial ecosystem. Our results showed that higher alpha diversity and lower network stability were observed in the proglacial ecosystem compared to the supra- and subglacial ecosystems. Furthermore, we found that the supraglacial ecosystem exported more high-abundance taxa, whereas the subglacial ecosystem released a more diverse taxa to the proglacial ecosystem.

Material and methods

Sampling

The Kuoqionggangri (hereafter ‘KQGR’, 29˚51′56″N, 90˚12′11″E) glacier catchment is located in the Mt. Nyainqen Tangglha in the south of the Tibetan Plateau (Fig. 1). Mean annual precipitation and temperature are 550 mm and −1.28 °C, respectively. Almost 90% of annual precipitation is concentrated in the summer from June to August44. There are seven small glaciers with areas ranging from 0.05 to 1.56 km2 in the catchment (Fig. 1), and at least nine glacier-fed lakes expanded on account of glacier meltwater and calving glaciers according to remote-sensing images. Diverse habitats exist in the terminus and foreland of the KQGR No. 1 glacier, such as snowpack, supraglacial streams and cryoconite, and proglacial lakes. Moreover, the temporary cave formed by glacier calving provides a window for investigating the subglacial ecosystem. Given the diverse habitats of the KQGR No. 1 glacier, we collected 42 samples from seven representative habitats of three glacial ecosystems (i.e., supraglacial, subglacial, and proglacial ecosystems) to reveal the similarities and differences in the microbial communities present within these ecosystems.

Fig. 1: Map and scheme of sampling sites.
figure 1

a Position of Kuoqionggangri catchment in the Lhasa basin in the Tibetan Plateau (source: Landsat imagery). b Schematic of glacier interconnected with streams and lake (source: Randolph Glacier Inventory version 6.0). c Map and photos of sampling sites. Sampling sites were labeled in remote-sensing image (source: Esri, Maxar, Earthstar Geographics) and sampling environments of each glacial ecosystem were displayed by filed photos. Sample abbreviation as fellows: SupC supraglacial cryoconite, SupS supraglacial stream, SubM subglacial meltwater, SubS subglacial sediment, ProW proglacial lake water, ProS proglacial lake sediment. Maps were created in ArcMap 10.2.

Full size image

All 42 samples were collected at KQGR No. 1 glacier (Figs. 1 and S1) and proglacial lake (4534 m2 area) between 5518 and 5627 m above sea level on 12th and 13th July 2020. All samples were collected from seven habitats of three glacial ecosystems, including snow (n = 8), cryoconite (n = 5, hereafter ‘SupC’), and stream (n = 7, SupS) in supraglacial ecosystem; sediment (n = 3, SubS) and meltwater (n = 3, SubM) in subglacial ecosystem; and sediment (n = 8, ProS) and water (n = 8, ProW) in the proglacial lake.

Snow samples were collected within 12 h after the snowfall event to avoid the development of cryoseston, and the upper 2 cm surface snow was removed for decontamination, and then placed in sterile WhirlPak (Nasco, WI, USA) bags by pre-cleaned steel scoop. Eight snow samples were collected at two sites with 2 m intervals in the glacier terminus. Five cryoconite granular samples were collected near the snow sample sites and placed at 50 mL sterile polyethylene centrifuge tubes using the pre-cleaned steel scoop. Supraglacial stream samples were collected using 500 mL sterile Nalgene bottles (Nalgene, NY, USA) in three sites about 30 m apart from upstream to downstream. We collected two samples at each site, with an additional sample collected at the third site.

Subglacial samples were collected in the cave formed by the glacier terminus calving event, which represents the recently exposed subglacial environment (Fig. S1). Owing to airflow in the subglacial cave, the exchange of heat between ice and air promotes the melting of the ice and the production of subglacial meltwater45. Meltwater generated on glacier surfaces also enters the subglacial drainage system. Subglacial meltwater samples were collected in three sites with 2 m intervals using 500 mL sterile Nalgene bottles (Nalgene, NY, USA). The subglacial sediments were covered by meltwater with an average depth of 3 cm, containing gravel, rock fragments, and deposited debris. After collected meltwater samples, subglacial sediment samples were collected into 50 mL sterile polyethylene centrifuge tubes using a cleaned scoop.

Proglacial lake water and sediment samples were collected at seven sites around the lake shore and one in the center of the lake. The lake water was collected using 500 mL sterile Nalgene bottles with ~20 cm depth and sediments were collected by the clean scoop in 50 mL sterile polyethylene centrifuge tubes. In particular, sediment samples in the center of the lake were collected using a grab sediment sampler. Clean and disposable nitrile gloves were worn during all sampling courses to avoid contamination. All samples were kept in an icebox during transport from the glacier to the laboratory within four hours of sampling, where samples were stored at −20 °C.

Nucleic acids extraction and sequencing

In the laboratory, each 1 L of thawed water of snow was filtered through a 47 mm polycarbonate membrane (0.22 μm, Millipore, USA) by a vacuum pump (Model GM-1.0 A, Tianjin Jinteng, China, 60 L per min, 0.08 MPa). Filters were cut into small pieces, 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 two hours before DNA extraction. DNA extraction of snow and water samples was performed by FastDNA® Spin kit (MP Biomedicals, Santa Ana, CA) following the manufacturer’s protocol. The frozen cryoconite and sediments were thawed at room temperature and 0.5 g of each sample was processed for DNA extraction using the FastDNA® Spin kit (MP Biomedicals, Santa Ana, CA) following the manufacturer’s protocol.

The hypervariable V4-V5 region of the 16S rRNA of bacteria was amplified by using primer pair 515 F (5′-GTGCCAGCMGCCGCGGTAA-3′) and 907 R (5′-CCGTCAATTCMTTTRAGTTT-3′). Template DNA from each sample was amplified in triplicate in a 50 μL reaction system with the following conditions: 30 cycles of denaturation at 94 °C for 30 s, annealing at 55 °C for 30 s, and extension at 72 °C for 30 s, with a final extension at 72 °C for 10 min. The amplified DNA was sequenced using the Illumina MiSeq platform (Illumina, San Diego, CA) with 250-bp paired-end reads at the Guangdong Magigene Biotechnology Co. Ltd. (Guangzhou, China). Reads generated in this study have been uploaded to the NCBI SRA database (BioProject accession number PRJNA1081793).

Bioinformatic analysis

Microbial bioinformatics was performed using QIIME 2 (version 2020.11)46. Raw sequence data were demultiplexed and quality-filtered using the q2‐demux plugin followed by denoising with DADA247 (via q2‐dada2). All amplicon sequence variants (ASVs) were aligned with MAFFT48 (via q2‐alignment) and used to construct a phylogeny with FastTree49 (via q2‐phylogeny). Alpha‐diversity metrics (Chao1 and Shannon index), and beta-diversity metrics (Bray‐Curtis dissimilarity) were estimated using q2‐diversity after samples were rarefied to 40,858 sequences per sample. Taxonomy was assigned to ASVs using the q2‐feature‐classifier50 against the Sliver 138 NR reference database51

Statistical analysis

To ascertain differences in community structure between habitats, we then used unconstrained ordinations (principal coordinates analysis; PCoA) to evaluate variability between habitats and ecosystems using Bray-Curtis distances. The significance comparisons were tested using a permutational multivariate analysis of variance (PERMANOVA) with the adonis function in the vegan (v2.6-2) package52. The significance of alpha diversity and SourceTracker results among glacial habitats and ecosystems was determined in R software (v4.2.1)53 using the Kruskal.test function from the R base package, followed by the DunnTest function from the FSA package for multiple comparisons.

Quantification of the community assembly process

We used a quantification framework based on a phylogenetic bin-based null model analysis (iCAMP) to quantify the relative importance of homogeneous selection (HoS), heterogeneous selection (HeS), homogenizing dispersal (HD), dispersal limitation (DL), and drift (DR) to the assembly of glacial ecosystems54. Selection includes both homogeneous selection (i.e., selective pressure resulting from the same environmental conditions) and heterogeneous selection (i.e., selective pressure resulting from different environmental conditions)55. Dispersal includes homogenizing dispersal (i.e., high dispersal rates homogenize community composition, leading to low compositional turnover) and dispersal limitation (i.e. low dispersal rates limit species exchange between local communities, leading to high compositional turnover)55. The assembly processes above mentioned were assessed by the β-net-relatedness index (βNRI) and Bray-Curtis-based Raup-Crick (RCBray) using the package iCAMP (v1.5.12) in R software (v4.0.3)56. The heterogeneous and homogeneous selections were quantified by βNRI higher than 1.96 and βNRI lower than −1.96, respectively. The relative importance of dispersal limitation and homogenizing dispersal was estimated by |βNRI | ≤ 1.96 but RCBray > +0.95 and RCBray < −0.95, respectively. Ecological drift was calculated by |βNRI | ≤ 1.96 and |RCBray | ≤ 0.95.

Co-occurrence network and robustness analysis

Correlation analysis was evaluated using the Sparse Correlations for Compositional Data (SparCC) algorithm in FastSpar (v1.0.0)57. First, ASVs present in over ~10% of samples were used to calculate correlation in FastSpar with default parameters. The significance of correlations was estimated using 1000 permutations and bootstrap. Next, significant correlations (P-value < 0.01 and |r | > 0.5) between ASVs were used to generate a sub-network of three ecosystems using the igraph package in R software (v4.0.3)56 and exported as a GML format network file for visualization in Gephi software. The network for each ecosystem was generated by keeping ASVs present in specific ecosystems as described by Ma et al.58. Topological features of each network were calculated in the igraph package, including node number, edge number, degree, betweenness centrality, average path length, density, and modularity. Yuan et al.43 successfully demonstrated that microbial network stability increases under climate warming by using network robustness assessment. The robustness of a network is defined as the proportion of remaining species in the network after random or targeted node removal43. In this study, the robustness was estimated by following the procedure of Yuan et al. for predicting the stability of microbial networks. In addition, we also analyzed the robustness of the proglacial network after removing the nodes shared between the supraglacial or subglacial ecosystem and the proglacial ecosystem to indicate the influence of supraglacial and subglacial species on the stability of the proglacial network.

SourceTracker

SourceTracker, a Bayesian mixing model59, was used to identify the proportions of each “sink” microbial community that is derived from each “source” microbial community. In this study, to quantify the contribution of the microbiome in supraglacial and subglacial ecosystems to the proglacial ecosystem, we used microbial communities in habitats of supraglacial and subglacial ecosystems as potential sources (i.e., snow, supraglacial streams, cryoconite, subglacial meltwater, and subglacial sediment), and microbial communities in the proglacial lake water and sediment as the final sink. The analysis workflow was described in the raw example. R script with default settings and alpha value set to 0.001 for simplifying model iteration. In addition, to verify whether supraglacial and subglacial habitats export different taxa downstream, we split the ASV table into four taxonomic groups (i.e., Alphaproteobacteria, Bacteroidota, Gammaproteobacteria, and others) for separate SourceTracker analysis. Groups Alphaproteobacteria, Bacteroidota, and Gammaproteobacteria represent the high-abundance taxa while the Others group represents the low-abundance taxa.

Results

Community taxonomic composition, diversity, and assembly

The microbial community was dominated by Bacteria, with Archaea present at a relative abundance of less than 1%. Among the top nine phyla, the phyla Proteobacteria and Bacteroidota accounted for ~88% relative abundance of the community (Fig. 2a). The community taxonomic composition of different glacial ecosystems showed slight variation (for Proteobacteria, we compared at the class level). Supra-, sub- and proglacial ecosystems were all dominated by Gammaproteobacteria (40.7%), followed by Alphaproteobacteria (29.5%) and Bacteroidota (18%) (Fig. 2a). Besides similar dominant taxa, the supraglacial ecosystem had more Firmicutes (4.5%) than the other two ecosystems, and the proglacial ecosystem observed more Acidobacteriota (4.2%). At the genus level, the top 14 genera accounted for 71% of the total relative abundance of the microbial community. Among these genera, Sphingomonas (39%) and Polaromonas (32%) were predominant across the three glacial ecosystems. The supraglacial ecosystem had more Acinetobacter (21%), the subglacial ecosystem had more Spirosoma (9%) and the proglacial ecosystem had more Brevundimonas (8%). Notably, the subglacial ecosystem had the highest relative abundance of methylotrophs, including the genera Methylobacterium, Methylotenera, and Methylorosula, which are closely associated with the subglacial environment60, compared to the other ecosystems (Fig. S2).

Fig. 2: Taxonomic composition of the microbial community at the phylum level across three glacial ecosystems and seven habitats.
figure 2

a Taxonomic composition of the microbial community across three glacial ecosystems. b Taxonomic composition of the microbial community across seven habitats. The phylum Proteobacteria was divided into class levels. Phylum that sums of relative abundance below top 9 merged as Others. Sample abbreviation as fellows: SupC supraglacial cryoconite, SupS supraglacial stream, SubM subglacial meltwater, SubS subglacial sediment, ProW proglacial lake water, ProS proglacial lake sediment.

Full size image

Adjacent habitats in the supraglacial ecosystem showed distinct taxonomic compositions (Fig. 2b). The snow was dominated by the Gammaproteobacteria (70%) and Alphaproteobacteria (15.5%), followed by Bacteroidota (12.3%). The cryoconite was dominated by Gammaproteobacteria (34.9%), followed by Bacteroidota (32.7%) and Firmicutes (14.2%). The stream was dominated by Alphaproteobacteria (40.4%) and Gammaproteobacteria (40.3%). As for habitats in subglacial and proglacial ecosystems, the taxonomic compositions of sediment and water were different. The subglacial meltwater and the proglacial lake water had more Alphaproteobacteria (SubM = 39.9%; ProW = 48.2%) than Gammaproteobacteria (SubM = 30.8%; ProW = 21.4%). Whereas subglacial sediments and proglacial lake sediments had more Gammaproteobacteria (SubS = 41.6%; ProS = 38.7%) rather than Alphaproteobacteria (SubS = 30.1% and ProS = 20.2%). Several unique phyla were observed high relative abundances (>3%) within a single habitat. For instance, the supraglacial cryoconite had more Cyanobacteria (5.2%) than other habitats, the proglacial lake had more Deinococcota (4.2%) and sediment had more Acidobacteriota (8.8%). At the genus level, in the supraglacial stream, subglacial meltwater, and proglacial lake water, the genus Sphingomonas was highly abundant, with an average relative abundance of 23%. While the subglacial and proglacial sediments contained a greater proportion of unclassified genera (26% and 16%, respectively), suggesting that these environments harbor diverse novel taxa. Additionally, habitat-specific genera with relative abundances exceeding 15% were identified. For instance, the snow contained a higher abundance of Acinetobacter (50%), supraglacial cryoconite had more Flavobacterium (17%), subglacial meltwater was characterized by an abundance of Spirosoma (18%) and proglacial lake water exhibited greater levels of Brevundimonas (16%).

The diversity of microbial community increased from the supraglacial to the subglacial and then to the proglacial ecosystem (Fig. 3a, b). The Chao1 and Shannon indices of the proglacial ecosystem were significantly higher than those of subglacial and supraglacial ecosystems (Kruskal–Wallis test, p < 0.001), while there was no significant difference between the subglacial and supraglacial ecosystems. Among the seven habitats, the Chao1 and Shannon index were significantly lower in the snow and higher in the proglacial lake sediments than in other habitats (Kruskal–Wallis test, both p < 0.05). No significant difference in alpha diversity among SupS, SubM, and ProW. The principal coordinates analysis showed that Snow, SupC, SubS, and ProS clustered separately, while SupS, SubM, and ProW clustered together (Fig. 3c). PCoA1 and PCoA2 together explained 42.76% of the community variation (Fig. 3c). Measures of between-sample diversity (beta diversity) among seven habitats revealed that the highest beta diversity of microbial communities was in the ProS and the lowest in the SubM (Fig. S3). PERMANOVA test also revealed significant differences in community structure between paired habitats (Supplementary Table S1). Habitat can explain 69% of the variation in total community structure (adonis R2 = 0.69; Fig. 3c). SupC exhibited the greatest difference from other habitats (average adonis R2 = 0.62), while SupS showed the smallest difference (average adonis R2 = 0.47).

Fig. 3: Alpha and beta diversity of different habitats and ecosystems.
figure 3

a Shannon index for taxonomy diversity in seven habitats and three ecosystems. b Chao1 index for taxonomy richness in seven habitats and three ecosystems; Boxes followed by different letters differed significantly at P < 0.05 (Kruskal–Wallis test). c Principal Coordinate Analysis (PCoA) plot of Bray-Curtis distances among seven habitats. Point shapes indicate ecosystem categorization. Sample abbreviation as fellows: SupC supraglacial cryoconite, SupS supraglacial stream, SubM subglacial meltwater, SubS subglacial sediment, ProW proglacial lake water, ProS proglacial lake sediment, SUP supraglacial ecosystem, SUB subglacial ecosystem, PRO proglacial ecosystem.

Full size image

Community assembly analysis showed that deterministic processes (i.e., heterogeneous selection and homogeneous selection) dominated supraglacial and subglacial ecosystems, whereas stochastic processes (i.e., dispersal limitation, homogenizing dispersal, and drift) dominated the proglacial ecosystem (Fig. S4; Supplementary Table S2). Comparing the relative contribution of these processes at the ecosystem level, the drift process (23%) had the highest contribution in the community assembly of the supraglacial ecosystem, homogeneous selection (57%) had the highest contribution in the community assembly of the subglacial ecosystem, and dispersal limitation (39%) had the highest contribution in the proglacial ecosystem. Among different glacial habitats, these ecological processes also showed distinct contributions. Snow was dominated by the drift process (61%). The supraglacial cryoconite, stream, and subglacial meltwater were dominated by homogeneous selection process (67, 79, and 76%, respectively). Subglacial sediments and proglacial lake sediments were both dominated by the dispersal limitation process (41%).

Contributions of the supraglacial and subglacial to proglacial ecosystem community

Source tracking analysis showed that both supraglacial and subglacial ecosystems were potential diversity sources of microbial community in proglacial ecosystems (Fig. 4). Community diversity of proglacial lake water is sourced from supraglacial (~48%) and subglacial (~40%) ecosystems, while the contributions between supraglacial and subglacial showed no significant difference (Fig. S5). For proglacial lake sediment, the contribution of the subglacial ecosystem (42%) was significantly higher than that of the supraglacial ecosystem (Kruskal–Wallis test; p < 0.001) (Fig. S5). More specifically, supraglacial streams and subglacial meltwater contributed 43.48 and 25.99% of the diversity of proglacial lake water, respectively (Fig. 4a). As for lake sediments, subglacial habitats altogether contributed ~42% of the total diversity (Fig. 4b). Overall, the supraglacial ecosystem essentially contributed only to the diversity of proglacial lake water, whereas the subglacial ecosystem contributed to the diversity of both proglacial lake water and sediment.

Fig. 4: Results of the SourceTracker-based Bayesian mixing model.
figure 4

a and b using proglacial lake water and lake sediment as a sink and other habitats as sources; c and d showing SourceTracker results of different taxa groups. proglacial lake water and lake sediment as a sink and other habitats as sources. Error bars represent the standard error of the mean. Sample abbreviation as fellows: SupC supraglacial cryoconite, SupS supraglacial stream, SubM subglacial meltwater, SubS subglacial sediment, ProW proglacial lake water, ProS proglacial lake sediment.

Full size image

To explore the influence of taxonomic differences on the source-sink relationship, we separately analyzed four taxonomic groups (i.e., Alphaproteobacteria, Bacteroidota, Gammaproteobacteria, and others) (Fig. 4c, d). In the Alphaproteobacteria group, supraglacial streams contributed the most diversity (43.13%) to the proglacial lake. Subglacial sediment contributes the most diversity (6.13%) to proglacial lake sediment. In the Bacteroidota group, supraglacial streams contributed most of the diversity (48.75%) to the proglacial lake. Subglacial sediment contributes the most diversity (7.25%) to proglacial lake sediment. In the Gammaproteobacteria group, supraglacial streams contributed the most diversity (29.13%) to proglacial lake water. As for proglacial lake sediment, subglacial meltwater contributed the most diversity (5.13%). In others group, which consists of diverse taxa with lower relative abundance, subglacial meltwater contributed the most diversity (51.75%). For proglacial sediment, only subglacial sediment contributed 3% of diversity, and the contribution of other habitats was negligible. The subglacial ecosystem showed a significantly higher contribution than the supraglacial ecosystem to proglacial lake sediment in all groups (Kruskal–Wallis test; p < 0.001; Fig. S6). Only in high-abundance groups, the supraglacial ecosystem had significantly higher contributions to proglacial lake water than the subglacial ecosystem (Kruskal–Wallis test; p < 0.001; Fig. S6). Overall, the supraglacial ecosystem exported more high-abundance taxa (i.e., Alphaproteobacteria, Gammaproteobacteria, and Bacteroidota groups) to the proglacial lake, whereas, the subglacial ecosystem exported more low-abundance taxa (i.e., Others group).

Community stability of three glacial ecosystems

Although only 5.1% of ASVs were common across all three glacial ecosystems (Fig. S7), these ASVs accounted for 81% of the relative abundance of the community. Among the common ASVs, 75% were associated with dominant taxa (i.e., Alphaproteobacteria, Gammaproteobacteria, and Bacteroidota), followed by 7% classified as Actinobacteriota and 6% as Cyanobacteria. At the genus level, aside from 17% of these ASVs being unclassified, the majority were identified as Sphingomonas (7%), followed by Flavobacterium (5%) and Hymenobacter (5%). Microbial co-occurrence network showed that 12% of total ASVs were significantly correlated with each other (Fig. 5a), a proportion greater than that of the common ASVs. This finding indicates that microbial interactions exist not only between different ecosystems but also within each ecosystem. Topological attributes of the networks were various among the three ecosystems (Supplementary Table S3). The supraglacial network showed the highest edges, degree, and density than other networks. The subglacial network showed minimal number of nodes and edges but had the highest modularity than other networks. The proglacial network had the highest nodes number and betweenness centralization but had the lowest degree and density. These results indicated the microbial interaction decrease from the supraglacial to the subglacial and proglacial ecosystems. Additionally, node taxonomy and edge attributes were slightly different among the three networks. Over 60% of nodes in the three networks belonged to the Alphaproteobacteria, Gammaproteobacteria, and Bacteroidota (Fig. 5a). However, the proglacial network had more nodes (6.6%) of the phylum Acidobacteriota. Although the percentage of negative correlation edges was less than 20% in three networks, the number of negative correlation edges increased from the supraglacial (824) to the subglacial (852) to the proglacial network (1062).

Fig. 5: Microbial co-occurrence networks and robustness analysis of three ecosystems.
figure 5

a Co-occurrence network of supraglacial, subglacial, and proglacial ecosystems from left to right. The sizes of the nodes (ASVs) are proportional to the number of connections. The taxonomy of each node is shown in different colors in the phylum level (phylum Proteobacteria was divided into class levels). Orange edges mean a positive correlation and green edges depict negative correlations. b Robustness measured as the proportion of taxa remained with 50% of the taxa randomly removed from the SparCC network of three ecosystems. Each error bar corresponds to the standard deviation of 100 repetitions of the simulation. The blue dashed line shows the mean robustness of the subglacial network. Error bars represent the standard error of the mean. Comparisons (Kruskal–Wallis test) among three ecosystems are indicated by significance level: *** (P < 0.001). c Robustness measured as the proportion of taxa remained with ASVs presented in supraglacial and subglacial ecosystems removed from the network of proglacial ecosystem. Sample abbreviation as fellows: SUP supraglacial ecosystem, SUB subglacial ecosystem, PRO proglacial ecosystem.

Full size image

The degree of nodes was also changed from the supraglacial and subglacial network to the proglacial network (Fig. S8). Besides the 239 nodes found in all networks, there were 60 and 64 nodes in the proglacial network that were shared with the supraglacial and subglacial networks, respectively. Our results showed that 45.7% of node degrees increased from the supraglacial and subglacial to the proglacial network, while 30.6% of node degrees remained unchanged. In all 239 common nodes, their degree was significantly higher in the proglacial network than in the supraglacial network (Kruskal–Wallis test; p < 0.001). Moreover, for the 60 supraglacial nodes and 64 subglacial nodes that shared with the proglacial network, their degree was significantly higher in the proglacial network (Kruskal–Wallis test; p < 0.001). We further simulated community stability under species extinction conditions by calculating robustness. When 50% of the nodes in each network were randomly removed, the supraglacial network had the highest robustness whereas the proglacial network had the lowest robustness (Kruskal–Wallis test; p < 0.001; Fig. 5b). When supraglacial nodes were removed, the robustness of the proglacial network decreased faster than when subglacial nodes were removed (Fig. 5c).

Discussion

Glacier melting under global warming increases the hydrological connectivity of different glacial ecosystems and habitats7. Our results revealed that hydrological connectivity influences the composition, diversity, and assembly of communities in different glacial habitats and ecosystems (Figs. 2, 3, and S5). Habitats with strong hydrological connectivity, i.e., supraglacial stream, subglacial meltwater, and proglacial lake water, showed more similar community composition, diversity, and assembly process than other habitats. Meanwhile, snow, supraglacial cryoconite, subglacial sediments, and proglacial lake sediments showed more distinct community composition, diversity, and assembly processes. Similar dominant taxa in different glacial habitats and ecosystems indicated species migration by hydrological channel (Fig. 2). In supraglacial streams, subglacial meltwater, and proglacial lake water, the genus Sphingomonas was found to be highly abundant. This genus comprises obligate aerobic microorganisms that have been reported in a wide variety of environments, such as glacier-fed rivers61, polar soil62 as well as marine sediments63. Habitat-specific taxa indicated that environmental filtering results in specific distribution of species in certain habitat types (Fig. 2). For instance, more genera Acinetobacter were found in snow habitats, and more Cyanobacteria were found in supraglacial cryoconite than other habitats. These habitat-specific taxa are also found in other glacier-related habitats, such as Arctic snow64 and Antarctic cryoconite65. However, when glacier hydrological connectivity increases, these habitat-specific taxa could export to downstream and increase the diversity of downstream ecosystems66.

The alpha diversity of the microbial community increased from supraglacial to proglacial ecosystem, implying that potential supraglacial and subglacial species migrated into proglacial ecosystems by hydrological channel. Previous research has shown that microbial diversity increased in glacier-fed ecosystems by collecting microorganisms from various sources, such as glacier-fed rivers61,67, proglacial lakes31,35,68, and proglacial soils69,70. Besides species migration, glacier meltwater brings sediments and nutrients into proglacial lake, which altering nutrient condition and turbidity. Previous studies have demonstrated that turbid proglacial lakes have higher microbial diversity and more light-inhibited taxa (e.g., Nitrospira) due to the rapid attenuation of UV radiation by suspended particles31. In addition, our results revealed the influences of hydrological connectivity in community diversity among different habitats (Fig. 3a). Although significant difference in alpha diversity was found in the ecosystem level, the alpha diversity of supraglacial stream, subglacial meltwater, and proglacial lake were close in habitat level (Fig. 3a, b), suggesting that hydrological connectivity could increase community similarity in the different glacial habitats9. These habitats connected by glacier meltwater also showed similar community structures due to homogeneous selection (Figs. 3c and S4). However, for other habitats with weak hydrological connectivity (i.e., supraglacial cryoconite, subglacial and proglacial sediments)71, their microbial community structure was distinct (Fig. 3c). Microbes of these habitats were relatively difficult to travel a long distance because they were attached to the sediments and their dispersal might be readily limited by available wind or water flow72. Furthermore, weak movement ability causes the long-time isolation that shapes distinct communities. Community assembly analysis also revealed that higher percentage of dispersal limitation process in supraglacial cryoconite, subglacial, and proglacial sediments, which suggests low dispersal rates (Fig. S4). All these results indicated that species migration from supraglacial and subglacial increased the diversity but decreased the environmental selection process of the proglacial ecosystem

Proglacial ecosystems are a potential sink for supraglacial and subglacial ecosystems; their community diversity and stability could be influenced by taxa from these ecosystems. Source-sink systems formed by supraglacial and proglacial ecosystems have been deeply studied in recent years31,73,74,75,76, but the relative contribution of subglacial ecosystems in this system was little known. Previous studies have reported the important contributions of supraglacial ecosystem to the increase the diversity of downstream ecosystems36,68. While both subglacial and supraglacial ecosystems play significant roles in shaping proglacial diversity, our source tracking analysis further revealed that the subglacial ecosystem contribute greater diversity to the proglacial ecosystem compared to the supraglacial ecosystem (Figs. 5a, b, and S5). High-abundance taxa in the community, such as Gammaproteobacteria, Alphaproteobacteria, and Bacteriodota, were mainly exported to the proglacial lake by the supraglacial stream, while low-abundance taxa were exported to the proglacial lake by subglacial water (Fig. S6). For the proglacial lake sediment, its diversity mainly derives from the subglacial ecosystem (Figs. 4b and S5). Together these results provide important insights into the existence of microbial mobility in three glacial ecosystems, with the proglacial ecosystem as the final collector receiving more diversity from subglacial than the supraglacial ecosystem.

In addition, network analysis can reveal the complex inter-interactions of different species in microbial ecology studies77. In this study, topological attributes of networks revealed the variation of microbial interactions and community stability from the supraglacial to the proglacial ecosystem (Supplementary Table S3; Fig. 4a). The species number and interspecies negative correlations both increased from the supraglacial to the proglacial network. More positive correlations were found in the supraglacial network, with the node being comprised with prevalence of photoautotrophic microorganisms in the supraglacial microbial community. Increased nodes and interspecies negative correlations in the proglacial network compared with supraglacial and subglacial network, could be explained by niche competition when microorganisms’ immigration. Early-arriving species would preempt available nutrients and weaken the survival ability of late arrivers78. Besides the change in microbial interactions, community stability indicated by the robustness of the network was also different among the three ecosystems (Fig. 4b). The proglacial ecosystem showed the lowest stability, suggesting that inputs from the supraglacial and the subglacial ecosystem could perturb the microbial community of the proglacial ecosystem. The perturbations caused by upstream inputs, such as suspended particles, alter the environmental characteristics of the proglacial ecosystem, thereby affecting microbial interactions34,79. Additionally, when upstream species enter this ecosystem, their interactions with resident species can also influence community stability. This result was in line with a higher influence of stochastic processes in the proglacial ecosystem (Fig. S4), indicating that the strong dispersal limitation can reduce the likelihood of species association and lead to unstable coexisting communities80,81. In addition, the robustness of the proglacial network decreased faster when nodes from the supraglacial ecosystem were purposefully removed (Fig. 5c). This result demonstrates that supraglacial species have a greater positive influence on the robustness of the proglacial network than subglacial species, suggesting the higher contribution of supraglacial species in maintaining the stability of the proglacial ecosystem than subglacial species. Additionally, these immigrants had significantly higher node degrees in proglacial network (Fig. S8). The degree of an individual node indicates its level of interaction with other nodes. A high degree suggests a central hub role in the network82. The elevated degree of supraglacial and subglacial nodes in the proglacial network indicated that more favorable conditions for these immigrants in proglacial ecosystem, such as higher temperature and nutrients. Our results also support the possibility that non-hub taxa of the supraglacial and subglacial become the hub taxa of proglacial ecosystem (Fig. S8), further indicating the importance of upstream species in maintaining the community stability of downstream ecosystems. It is crucial to acknowledge that our subglacial samples may not accurately reflect the strict subglacial environment (i.e., dark and anaerobic) due to its exposure to aerobic and light conditions. However, the microbial community in the basal zone of the glacier terminus has often been used to represent or infer the characteristics of the subglacial ecosystem24,69,83, as microbial activity was still constrained by the harsh cold environment. Additionally, methylotrophs associated with the subglacial environment60 had higher relative abundance in the subglacial ecosystem than in other glacial ecosystems (Fig. S2), which supported the representativeness of our subglacial samples. Overall, our results support the hypothesis that the contributions of the supraglacial and subglacial to the proglacial ecosystem are different. Our results highlight that supraglacial and subglacial ecosystems respectively influence the stability and diversity of microbial community in the proglacial ecosystem, but further studies with more complete information about microbial metabolic potential and activity are needed to clarify the mechanisms. Importantly, all these results together provide key information to predict possible shifts in community diversity and stability with climate change in the future.

Conclusion

Rapid glaciers shrinking accelerates the microbial downstream migration. We demonstrate that microbial community diversity, taxonomic composition and stability were distinct among seven habitats in three glacial ecosystems. However, habitats with strong hydrological connectivity exhibited similar community structures, suggesting microbial downstream dispersal and indicating the integration of communities across different habitats. In addition, our results highlight that supraglacial and subglacial ecosystems have a different impact on the proglacial ecosystem, with the former exporting mainly high-abundance taxa and being more conducive to proglacial ecosystem stabilization, while the latter releasing a more diverse range of low-abundance taxa, increasing proglacial ecosystem diversity. Against the background of accelerated glacier retreat, more ancient or rare organisms from glaciers may be detected in proglacial ecosystems, which could help to improve our understanding of glacial microbial diversity and predict their ecological effects. Collectively, our findings contribute to improve the understanding of microbial shift and the connection between glacial and glacier-adjacent ecosystems and to predict the impacts of glacier retreat on biodiversity and ecosystem stability.