Predicting climate-change impacts on the global glacier-fed stream microbiome

Domain: earth_science
Published: 2025-02-01
Editorial outcome: minor_revision

Abstract
The shrinkage of glaciers and the vanishing of glacier-fed streams (GFSs) are emblematic of climate change. However, forecasts of how GFS microbiome structure and function will change under projected climate change scenarios are lacking. Combining 2,333 prokaryotic metagenome-assembled genomes with climatic, glaciological, and environmental data collected by the Vanishing Glaciers project from 164 GFSs draining Earth’s major mountain ranges, we here predict the future of the GFS microbiome until the end of the century under various climate change scenarios. Our model framework is rooted in a space-for-time substitution design and leverages statistical learning approaches. We predict that declining environmental selection promotes primary production in GFSs, stimulating both bacterial biomass and biodiversity. Concomitantly, predictions suggest that the phylogenetic structure of the GFS microbiome will change and entire bacterial clades are at risk. Furthermore, genomic projections reveal that microbiome functions will shift, with intensified solar energy acquisition pathways, heterotrophy and algal-bacterial interactions. Altogether, we project a ‘greener’ future of the world’s GFSs accompanied by a loss of clades that have adapted to environmental harshness, with consequences for ecosystem functioning.

Introduction
Predicting the impacts of climate change on biodiversity has become a mainstay in ecological research1,2. Despite the intricate relationships between microorganisms and climate, climate change microbiology is still in its infancy, which is particularly true for forecasting responses of entire microbiomes to climate change3. This is surprising, given the predictive power encoded in microbiomes, which integrate past and current environmental conditions and drive key ecosystem functions4. Despite this, relatively few studies have predicted climate impacts on microbiome structure and function—mostly in the ocean5,6 and soils7,8. Today, no such study exists for stream and river microbiomes, which contrasts their relevance for global biogeochemical cycling9 and ecosystem services10.

Glacier-fed streams (GFSs) initiate the flow of water for some of the world’s largest river systems and provide water resources to large human populations11, but are also most vulnerable to climate change. The shrinking of glaciers not only threatens water availability but also fundamentally alters the GFS environment, putting their ecological communities at risk. Recent studies on invertebrates suggest that species adapted to the harsh environmental conditions of GFSs become increasingly imperilled as glaciers disappear, and their distributions are reduced to cold-water refugia12,13. This biodiversity at higher trophic levels is to a large extent sustained by microbial biofilms that coat the GFS streambeds and regulate key ecosystem processes, including metabolism and nutrient cycling14. The genomic repertoire of biofilm-dwelling microorganisms allows them to cope with the highly selective GFS environment (e.g. near-freezing water temperatures, high UV radiation, and ultra-oligotrophy), and to seize opportunities when they become available at punctuated periods of the year15,16. Recent studies based on space-for-time substitution approaches suggest that climate change-induced glacier shrinkage may stimulate organic matter decomposition by microorganisms in GFSs17,18 and shift microbial energetics from chemolithoautotrophy to heterotrophy15,17. However, predictions of how the GFS microbiome structure and function will change under projected climate change scenarios are still lacking.

Here, we develop a hierarchical, machine-learning-based modelling framework to forecast how the structure and function of the global GFS microbiome will change in response to climate change under different greenhouse gas emissions scenarios (Shared Socioeconomic Pathways; SSP19) for the years 2070–2100. We leverage 2333 bacterial metagenome-assembled genomes (MAGs, strain-level resolution, 99% average nucleotide identity) and 6226 KEGG orthologs, alongside glaciological and environmental data, from GFSs sampled by the Vanishing Glaciers project across the Caucasus Mountains, the European Alps, the Andes, Himalayas, Pamir and Tien Shan, Rwenzori Mountains, Scandinavian Mountains, New Zealand Southern Alps, and Southwest Greenland (Fig. 1A). In each of the 164 GFSs, we captured the deglaciation history since the end of the Little Ice Age period, mirroring our forecasting horizon (Supplementary Fig. 1). Our forecasts show how the global GFS microbiome, including its biomass and biodiversity, as well as structure and function will shift under future climate scenarios.

Results
Climate-induced changes in the glacier-fed stream environment
The structure and function of stream biofilms are shaped by a suite of environmental factors (e.g. flow-induced hydraulics, temperature, resources), which are changing predictably as glaciers shrink14. In fact, glacier shrinkage alters flow and temperature regimes in GFSs, as well as sediment loads, with implications for turbidity and resources (e.g. phosphorus, nitrogen, organic carbon)17,20–22. Given the influence the environmental template has over the GFS microbiome23, we first modelled the future GFS environment (2070–2100), considering SSP climate change scenarios, of which we here report the results for SSP3 (results for SSP1 and SSP5 are provided in Supplementary Information). Our modelling framework is rooted in glaciological and climatological projections based on the Global Glacier Evolution Model (GloGEM)24 and Climatologies for Earth’s Land Surface Areas (CHELSA)25 at high spatial resolution, respectively. Combining these projections with measured data of the GFS environment along the chronosequence (i.e. leveraging differences between up- and downstream reaches for each GFS) (“Methods”, Supplementary Fig. 1), we first predict how the GFS environment will change in response to glacier retreat by the end of the century.

Our predictions suggest that median streamwater temperature and electrical conductivity (a proxy for ion concentration) will increase by 306.7% (IQR: 87.9–633.1%) and 88.2% (IQR: 51.6–130.7%), respectively, while median streamwater turbidity will decrease by 44.4% (IQR: 31.6–71.7%) across all study GFSs (Supplementary Fig. 2). Median concentrations of soluble reactive phosphorus and dissolved inorganic nitrogen will decrease by 14.1% (IQR: 9.6–27.3%) and 11.5% (IQR: 5.2–16.3%), respectively, and pH will be lowered by 2.8% (IQR: 1.6–4.2%) (Supplementary Fig. 2). Response curve analyses suggest that the magnitude of changes depend on glacier area (Supplementary Fig. 3). Consistent with previous studies and conceptual models20,26, these findings suggest reduced erosion and weathering capacity of shrinking glaciers, leading to decreases in streamwater fine suspended sediments and soluble reactive phosphorus concentrations. Changes in these selective constraints are likely to have significant impacts on the microbiome.

Biodiversity shifts with the greening of glacier-fed streams
Benthic microbial biomass is key to stream ecosystem functioning, as it fuels the food web and regulates ecosystem energetics and nutrient cycling14. Our projections suggest significant increases in benthic chlorophyll a (339.7%; IQR: 183–852.2%), a proxy for algal biomass, and bacterial abundance (88.5%; IQR: 60.4–150.2%) until the end of the century (Fig. 1B, C). This projection corroborates recent evidence suggesting that GFSs become ‘greener’ because of the fading capacity of glaciers to generate fine sediments, thereby reducing light limitation in GFSs17. In fact, fine suspended sediments render GFSs turbid, which reduces light availability for photosynthesis and increases physical abrasion, a notion that is supported by a negative correlation between projected values of streamwater turbidity and benthic chlorophyll a (rho = −0.96, p < 0.001). Yet, the projected benthic biomass of GFSs is still low compared to high-alpine streams without glacial influence27, which points to the persistence of constraints in GFSs other than those related to flow and turbidity.

Despite the harsh GFS environment, resident biofilms host a distinct and diverse microbiome15,23,28. Our projections show that this diversity (expressed as Shannon H’) will increase by 6.2% (IQR: 4.7–8.9%) under SSP3 (Fig. 1D). Positive relationships of both, bacterial Shannon diversity (rho = 0.69, p < 0.01) and bacterial abundance (rho = 0.98, p < 0.001), with chlorophyll a supports the notion that increased primary production sustains higher microbial biodiversity. Underpinning mechanisms may include elevated environmental stability and energy availability, the latter being a well-known driver of biodiversity29. Hence, we predict that future increases in algal biomass, and thus resource availability, will impact GFS microbiome diversity.

Both deterministic and stochastic assembly processes imprint upon the phylogenetic structure of microbial communities30. Environmental filtering, for instance, can lead to phylogenetic clustering. In GFSs, deterministic environmental selection favours microdiverse clades, which greatly contribute to microbiome diversity31. We posit that weaker selective constraints under future climates will alter microbial community assembly processes. To explore this emerging property, we modelled the mean nearest taxon distance and mean phylogenetic distance, reflecting relatedness at shallow and deep phylogenetic branching, respectively. We find median values of mean nearest taxon distance to significantly (Generalized Additive Model, GAM p < 0.001) increase by 3.5% (IQR: 2.2–5.8%, Fig. 1E), suggesting terminal phylogenetic clustering will diminish in the future GFS microbiome. Correlation between predicted mean nearest taxon distance and chlorophyll a (rho = 0.97, p < 0.001) suggests increased resource availability may contribute to the reduction of selective constraints. While the ability to use different carbon substrates can promote divergence within microdiverse clades32, microdiversity in GFSs is concentrated in clades utilising chemolithotrophic energy pathways31. Reduced environmental selection may erode microdiversity in GFSs, with yet unknown consequences for the stability and resilience of the GFS microbiome33. Similarly, the mean phylogenetic distance is predicted to increase by 3.2% (IQR: 1.7–6.3%), pointing towards deeper-branching changes to the future GFS microbiome phylogeny. In line with the predicted increase in alpha diversity, we attribute such a deeper-branching expansion of the phylogeny to the establishment of novel taxa and lineages in future GFSs.

Climate change shifts strain distributions
Species distribution models are commonly used to forecast climate-change impacts on the structure and diversity of ecological communities34. To predict the impacts of climate change, and integrate glaciological and environmental controls on the GFS microbiome, we built individual models for each of the 2333 strain-level MAGs in a strain distribution model framework (“Methods”, Supplementary Fig. 1). For this, we used a combination of climatic, glaciological, and mineralogical data, as well as forecasts of streamwater physico-chemistry. Overall, these models predicted strain abundance with satisfactory accuracy (median cross-validation R2 = 0.25; IQR: 0.13–0.36; Supplementary Fig. 4).

We first assessed the importance of covariates selected as predictors to identify the main drivers of strain abundances (Supplementary Fig. 5). Across all strains, streamwater electrical conductivity, pH, and temperature were the most important predictors, along with latitude. Distance to the glacier, annual snow cover, and bioclimatic variables (reduced by Principal Component Analysis) were also identified as important predictors (Supplementary Fig. 5). Interestingly, chlorophyll a did not rank among the strongest predictors of strain abundances. Together, these predictors highlight the direct controls of glacial meltwaters and underlying bedrock geology on the abundance of microorganisms in GFSs. Moreover, phylogenetically closely related strains shared similar predictors compared to less related strains (Spearman correlation, rho = −0.15, p < 0.0001, Supplementary Fig. 6). This apparent niche conservatism supports our modelling approach.

Next, using forecasts of these predictors, we assessed future abundance distributions for each strain. The majority (i.e. 64.7%) of the 2333 strains are expected to increase in abundance, with only 5.3% remaining unchanged (Fig. 2A and Supplementary Table 3). This overall gain in abundance aligns with our independent forecast of increasing bacterial abundance. Overall, Gamma- and Alphaproteobacteria, which numerically dominate the present-day GFS microbiome15, are projected to experience the largest increase in abundance (Fig. 2C). This increase in abundance can, at least partially, be attributed to the large number of strains in these classes. However, strains that currently occur at low abundance (i.e. lower half of the abundance distribution) are projected to increase disproportionately more in abundance compared to strains that are currently abundant in GFS microbiomes (Fig. 2B). Among them, the class of Gemmatimonadetes, which is known to form close associations with freshwater algae35 and Paceibacteria, known for their parasitic or symbiotic lifestyles36, revealed the largest relative increase in abundance. We hypothesise that Paceibacteria may be promoted by a larger size of the host pool, putatively related to elevated microbial biomass in future GFSs.

We also found that 30% of the strains will decrease in abundance (Fig. 2B), further pointing to the reorganisation of the GFS microbiome. Unlike covariates associated with strains predicted to increase in abundance, annual snow cover duration, distance to the glacier, cryospheric influence in GFS are facing a reduction in abundance in the future. This is underlined by the fact that changes in abundance were more similar for phylogenetically closely related strains than for distantly related strains (i.e. phylogenetic signal in log2-fold abundance change: lambda = 0.88, p < 0.001, Fig. 2A and Supplementary Table 4). The finding that phylogenetically related strains decrease in abundance reflects their shared evolutionary history and suggests that fundamental ecological niches in GFSs will undergo major transformations as glaciers shrink. This raises concerns, as it may mean that climate change imperils entire clades in GFSs rather than individual strains. This concern is substantiated by the fact that 26.6% of the strains projected to decrease in abundance were members of monophyletic clades, of which all members are projected to decrease in abundance (Supplementary Table 5). Notably, some of the largest clades (e.g. Ferruginibacter, Lacisediminihabitans, Acetobacteraceae) are hallmarks of cryospheric ecosystems37–39. Taken together, these findings suggest that the microbiome of the world’s GFSs will experience a profound compositional and phylogenetic restructuring with climate change. Although overall diversity of the GFS microbiome may expand due to the availability of novel niches associated with increased algal biomass the loss of microdiversity in putatively well-adapted clades likely has important ramifications for community stability and resilience.

Genomic traits associated with decreased abundance
To gain more mechanistic insights into the projected trajectories of strains, we next linked genomic traits to predicted changes in abundance. We found that strains predicted to decrease in abundance have generally smaller genomes (GAM, fixed effect estimate = −0.2 mbp, p < 0.01) but encode more KEGG (Kyoto Encyclopaedia of Genes and Genomes) ortholog groups (KOs; GAM, mean difference = 4.02, p < 0.001), compared to strains that will increase or stay invariant in abundance (Fig. 3A). Elevated metabolic diversity (i.e. higher KO numbers), despite reduced genome size, points towards genome optimisation, which can be explained by reduced KO redundancy. Indeed, estimating KO redundancy as the ratio of unique KOs to all KOs on genomes, we found that strains predicted to decrease in abundance have reduced KO redundancy (GAM, mean difference = −0.0136, p < 0.001).

Together, this suggests that despite a broad metabolic repertoire, smaller genomes are favoured under current climatic and environmental conditions in GFSs, putatively reflecting the adaptation to ultra-oligotrophy and the unstable GFS environment. This notion is supported by the observation of a large degree of mixotrophy among GFS bacteria, supposedly enabling them to exploit varying energy sources (e.g. solar radiation and organic carbon15). However, these traits may no longer be favoured in GFSs under a future climate.

Given the predicted decline in abundance of entire clades, we next focused on KOs associated with strains that are projected to decrease in abundance. To this end, we constructed random forest classifiers and used a leave-one-out approach accounting for phylogenetic structure (“Methods”), identifying 408 KOs that characterise these strains (p < 0.05, feature importance quantile > 0.95). KOs associated with biofilm formation and cold adaptation were particularly prevalent (Supplementary Table 6). These include, for instance, the sec-independent protein translocase (TatB) which has previously been linked to cold-shock in Shewanella oneidensis40, and DNA topoisomerase I (topA), whose mutants become cold sensitive in Escherichia coli41.

Using enrichment analysis, we found KEGG categories related to carbohydrate metabolism (Fisher test, OR = 2.20, p < 0.001) and energy metabolism (Fisher test, OR = 1.80, p < 0.001) to be significantly enriched among strains predicted to decrease in abundance (Fig. 3B and Supplementary Table 7). The prevalence of these KEGG categories may be linked to changing glacial and terrestrial carbon sources within catchments42,43. In fact, as GFSs become ‘greener’, carbon sources may shift towards more predictable algal-derived organic carbon17. Consequently, MAGs encoding chemolithoautotrophic metabolisms may become less prevalent whereas metabolisms such as aerobic chemoorganoheterotrophy may become more prevalent15. These results thus provide evidence for an alleviation of the multi-faceted selective constraints on future GFS microbiomes, ultimately altering their function.

Alteration of the functional potential
The predicted compositional reorganisation of the GFS microbiome is expected to affect its metabolic repertoire. Besides resolving genomic traits and functions associated with strains that will decrease or increase in abundance in the future, we modelled shifts in microbiome functional potential by building 6226 individual models of KO abundances. Compared to the strain-resolved models, these models yielded generally lower predictability (median cross validation R2 = 0.01, IQR: 0.0009–0.053), which is intuitive because of the invariance of core functions. However, 45.7% of KOs are predicted to increase, whereas 38.4% were predicted to decrease in abundance (Fig. 4). Pathways with the largest proportions of KOs forecasted to increase in abundance included photosynthesis (88% of KOs in this pathway are predicted to increase) and photosynthesis antennae (88%). We posit that this putatively reflects the predicted increase in light availability because of reduced turbidity in future GFSs. Several more specific metabolic pathways are also expected to increase in abundance, including selenocompounds metabolism (61.5%), alanine, aspartate and glutamate metabolism (62.7%), taurine and hypotaurine metabolism (60%), and nitrogen metabolism (58.7%) (Fig. 4). Pathways enriched in KOs predicted to decrease in abundance, and thus expected to play a reduced role in the future GFS microbiome, were related to protein processing (62.1%), sphingolipid metabolism (61.9%), replication and repair (60.4%), thiamine metabolism (59.1%), peptidoglycan biosynthesis (56.7%), fatty acid degradation (54.5%), and cell growth (51.7%) (Fig. 4). Taken together, the eminence of microbiome functional turnover suggests a substantial functional reorganisation of the future GFS microbiome with potential impacts on ecosystem functioning (e.g. decomposition and nutrient cycling).

In conclusion, we provide global-scale projections of the GFS microbiome for climate change scenarios. Environmental constraints that vary little across GFSs worldwide, but impose strong selective forces on the microbiome, allowed us to forecast changes in the environment and linked to this, of microbiome structure and function. Overall, we found that climate-induced glacier shrinkage relieves GFS microbiomes from environmental constraints, putting at risk entire clades of microorganisms that have adapted to the extreme environment. Our projections also reveal that microbiome functions will shift, particularly towards light acquisition pathways and translation, which is in line with the ‘greening’ of GFSs. We acknowledge potential caveats inherently related to species distribution models, including potentially missing variables and overfitting. We argue that a priori knowledge of the GFS environment and its microbiome, as well as cross-validation and ensemble modelling collectively make our projection architecture robust. Ultimately, future expeditions and long-term observations will allow a better understanding of climate change impacts on the microbiome of a rapidly vanishing ecosystem.

Methods
Study sites and sample collection
GFSs were sampled by the Vanishing Glaciers project between January 2019 and July 2022 (Supplementary Fig. 1 and Supplementary Table 1), and the resulting stream parameters’ dataset is also described in ref. 15 was used for this study (physicochemical parameters, stream water characteristics, benthic chlorophyll-a and bacterial abundance). Our global sampling included the European Alps, Scandinavian Mountains, Himalayas, Pamirs and Tian Shan, Ecuadorian and Chilean Andes, Southwest Greenland, Alaska Range, Caucasus, Rwenzori Mountains, and the New Zealand Southern Alps. For the sake of comparability, GFSs were predominantly sampled in spring or autumn during ‘windows of opportunities’16, which are characterised as periods either directly preceding or following peak glacier melt. For each GFS, two reaches were sampled: one as close as possible to the glacier terminus (median distance to glacier snout = 76 m, IQR = 29–301), and one located further downstream (median distance = 706 m, IQR = 336–1280). This approach allowed us to capture changes in glacier influence over the two stream reaches, creating a space-for-time substitution design. Space-for-time substitutions are well suited for communities shaped by environmental selection, but require communities to be part of the same ecological trajectory. This constrained the location of downstream reaches to be situated close to the terminal moraine, in many mountain regions hallmarks of the last glacial maximum (Little Ice Age). Because of the poorly developed soils and lack of vegetation in catchments above the terminal moraine, the influence of glacial meltwater (e.g. low temperature, turbidity, inorganic nutrient and organic carbon biogeochemistry) determines the environmental template in these reaches. At both, up- and downstream reach, we sampled three separate sediment patches (<10 m distance) with flame-sterilised devices. At each patch, we collected sandy (250 μm–3.15 mm; Retsch) sediments from the upper 5 cm of the benthic zone. Sediment samples were transferred to sterile cryovials, immediately flash-frozen in liquid nitrogen, and subsequently stored at −80°C prior to and following shipping to Switzerland for DNA extraction and biomass analyses.

Streamwater and sediment physicochemical parameters
We measured streamwater temperature, pH, specific conductivity, and turbidity using a portable MultiLine Multi 3630 IDS (WTW) metre. Streamwater samples for nutrient analyses were filtered through pre-combusted GF/F filters (Whatman) into acid-washed Nalgene HDPE bottles and frozen within 48 h. Samples for the determination of DOC concentration were filtered identically, but into pre-combusted glass vials and kept at 4°C. A LaChat QuikChem 8500 flow injection analyser was used to measure concentrations of ammonium (N-NH4+; QuikChem method 10-107-06-3-D), nitrate (N-NO3−, QuikChem method 10-107-05-1-C), and soluble reactive phosphorus (P-SRP; Method 10-115-01-1-M). Inorganic forms of nitrogen were combined to quantify dissolved inorganic nitrogen (DIN). A Sievers M9 TOC Analyser (GE) was used to measure DOC concentrations. Major cation and anion concentrations were measured from sterile-filtered (0.2 μm, Millipore) streamwater samples (stored in the dark at 4°C) using a Metrohm 930 Compact IC flex system. Benthic sediment mineralogy was determined using an X-TRA Thermo ARL Diffractometer (errors varied between 5% and 10% for the phyllosilicates and 5% for grain minerals). Relative abundances of the main mineral groups (clays, quartz, feldspar, and calcite) were computed from counts of mica, chlorite, amphibole, feldspars, calcite, and quartz.

Benthic chlorophyll a and bacterial abundance
We measured benthic chlorophyll a concentrations as a proxy for benthic algal biomass. For this, we added 90% ethanol to sediment samples, and incubated the samples in a hot water bath (78°C) for 10 min, followed by incubation in the dark at 4°C for 24 h. Then, samples were vortexed, and centrifuged and the supernatant was read in a high-sensitivity plate reader (BioTek Synergy H1) at 436/680 ex/em wavelengths. Chlorophyll a concentrations were calculated from standard curves obtained using a spinach standard (Sigma Aldrich) and normalised to a dry mass of sediment. Bacterial abundance was determined using flow cytometry as described previously31. Briefly, 2.5–3 g of sediment were fixed using a paraformaldehyde/glutaraldehyde mixture and kept frozen. In the laboratory, sodium pyrophosphate (final concentration of 0.025 mM) was amended and microbial cells were detached from sediments using vigorous shaking and sonication (Sonifier 450, Branson, 1 min, 60% duty cycle, output 5). Extracts were briefly spun, the supernatant 10-fold diluted and stained using SYBR Green and analysed on a NovoCyte flow cytometer (ACEA Biosciences) equipped with a 488 nm laser. Results were corrected for dilutions and normalised to sediment dry mass.

Climatic data
Climatology data at high spatial resolution was obtained from the CHELSA database (version 2.1)25. CHELSA provides both climatic and bioclimatic data typically used in species distribution modelling approaches at high spatial resolution (~1 km) based on a global downscaling approach. The database is based on an ERA-Interim climatic reanalysis and contains future projections of climatology for Shared Socioeconomic Pathways (SSP)19. The data was obtained using Python and the database’s API, and processed with the rasterio (v1.3.8) and gdal (v3.7.0) Python libraries44,45. GPS coordinates of all sampling locations were used to identify the corresponding grid cells of the (bio)climatic dataset. The specific sampling months were used for monthly parameters. For future projections, five different institutional models were downloaded and combined by averaging (gfdl-esm4, ukesm1-0-ll, mpi-esm1-2-hr, ipsl-cm6a-lr, and mri-esm2-0). Because values in the database were only available up to 2010, linear extrapolation based on the 1981–2010 period (per month) was performed to obtain data matching the 2019–2022 sampling period (Supplementary Fig. 7). In addition, climatic data for the 2070–2100 time period for three SSPs scenarios corresponding to different greenhouse gas emission trajectories were obtained. These scenarios included SSP1, SSP3, and SSP5, encompassing a range of potential future climate outcomes. We report median change and interquartile ranges for the site-specific changes for all parameters. To test for significant shifts, we conducted Wilcoxon tests of the future projections minus the present conditions. All statistical analyses were performed using R (v4.3.1)46 and the tidyverse R package suite (v2.0.0)47. Figures were prepared using the ggplot2 (v3.4.2)48 and ggridges R packages (v0.5.4)49.

Glaciological data
The future evolution of all mountain glaciers is available through the Global Glacier Evolution Model (GloGEM)24. The model is initialised with present-day glacier extents and computes changes in snow accumulation and melt, as well as changes in glacier length based on an ensemble of Global Circulation Models using different greenhouse-gas emission scenarios. The model has been calibrated to match observed mass changes at the scale of every individual glacier globally50. For this study, we extracted the distance of sampling locations to the glacier terminus, as well as the area of the glacierized surface for the time period 2070–2100. We report median change and interquartile ranges for the site-specific changes for all parameters. To test for a significant shift, we conducted Wilcoxon tests of the future projections minus present-day observations.

Metagenomics
Metagenomes were sequenced for 155 sediment samples covering all visited mountain ranges except the Alaska Range (these samples were collected last by the Vanishing Glacier Project and were not included in metagenomic sequencing runs). DNA extraction, purification, library preparation, sequencing and metagenome assembly steps were performed as described in ref. 51. Briefly, 5 g of GFS sediments were treated using an optimised phenol:chloroform-based extraction method subsequently followed by an ethanol precipitation step. This protocol yielded on average 50 ng of DNA per sample which was used for library preparation using the NEBNext Ultra II FS library kit, which included 6 PCR cycles. Sequencing library quality was screened using Qubit (Invitrogen) and Bioanalyzer (Agilent). Subsequently, libraries of app. Sixty samples were pooled and sequenced at the Functional Genomics Centre Zurich on a NovaSeq (Illumina, 150 bp) using S4 flow cells.

The metagenomic sequence data generated for ref. 15. was used, and was processed using the Integrated Meta-omic Pipeline workflow (version 3.0; commit# 9672c874)52. Briefly, adaptor trimming from reads is followed by an iterative assembly using MEGAHIT53 and metaFlye54. Reads from each individual sample were assembled into contigs and subjected to multi-coverage binning to obtain MAGs. For each individual assembly, we mapped the reads of the 5 geographically closest samples (Euclidean distances of gps coordinates) using BWA-mem (v0.7.17)55. To reduce computation time, we removed sequences in the assembly of <1.5 kbp. Subsequently, 10% of the pre-processed reads were randomly selected before mapping with seqtk (v1.3)56. We then used MetaBAT2 (v2.15)57, CONCOCT (v1.1.0)58 and MetaBinner (v1.4.3)59 using default parameters to obtain bins. DASTool (v1.1.4)60 was then employed to generate a non-redundant set of bins using a score threshold of 0.3. The quality of bins was assessed using CheckM2 (v1.0.1)61. Bins from all samples (including the ones generated by IMP3) with completeness of more than 50% (n = 12,599 bins) were selected for further analyses. We used MDMCleaner (v0.8.3)62 to reduce contamination. Finally, after rerunning CheckM2 to get final estimates of completeness and contamination, we used dRep (v3.2.2)63 to dereplicate bins using minimum completeness of 70% and maximum contamination of 10% and an ANI of 99%. This led to 2868 MAGs (strain-level dereplication). Functional annotation of the MAGs was performed with eggNOG-mapper (v2.1.9)64 against eggNOG v5.0 after obtaining coding regions with prodigal (v2.6.3)65. The coverage of the MAGs was estimated by mapping the reads of the samples to the genomic contigs using CoverM (v0.6.1, available at https://github.com/wwood/CoverM) using the trimmed_mean parameter. We normalized the coverage by similarly mapping the reads on the recA gene (K03553). After removing low-abundance strains (prevalence lower than 20% at a 10× recA coverage threshold), 2333 were selected for the strain distribution modelling. GTDB-Tk (v2.1)66 was used to determine the taxonomy of the MAGs. We used the concatenated alignment of 120 ubiquitous single-copy proteins created by GTDB-Tk to de novo generate a phylogenetic tree using FastTree2 under the WAG model of protein evolution with gamma-distributed rate heterogeneity.

Modelling environmental parameters and microbial biomass
Models of streamwater temperature, turbidity, electrical conductivity, pH, soluble reactive phosphorus (SRP), and dissolved inorganic nitrogen (DIN) were built using climatic, bioclimatic, glaciological and geological parameters as covariates, chosen by automated feature selection (Supplementary Fig. 3). The same modelling approach was used to model biomass parameters (chlorophyll-α, bacterial abundance) and diversity metrics (Shannon H’, mean nearest taxon distance, mean phylogenetic distance). Response variables and covariates were log-transformed where necessary to improve the normality of residuals, adding a constant equal to half of the smallest non-zero value. Before model fitting, all variables were z-transformed (i.e. scaled to mean and units standard deviation). All generalised additive models (GAM) models were created using the Gaussian family function and the bam function of the mgcv package (v1.9_0)67. GAMs have been successfully used to model and predict future changes in environmental sciences and ecology, especially for spatially structured data7,68,69, and have also been applied to stream ecosystems70,71. Here, GAMs were used to model stream parameters and biomass in order to account for spatial autocorrelation at the regional scale using a spatial spline (formula: s(latitude, longitude, bs=‘sos’, m=1, k=−1)). Feature selection was performed by building individual GAM models with a spline for each covariate (k=3, bs=‘ts’) along with the spatial spline and were then ranked by −log(p-value) of the covariate spline. The top three variables were retained for a final model. This procedure was repeated on each of the nine models for each cross-validation fold. A final model was then built with the spatial spline and a spline for the three top variables with the following parameters: k=3, bs=‘ts’. These parameters allowed for non-linear relationships while the small number of knots and the penalisation on the spline were added to control for smoothness and to avoid overfitting.

The models were validated using 10-fold cross-validation, meaning that several times the models were trained on 90% of the glacier-fed streams, and tested on the remaining data. The performance of the models (measured as the cross-validation R2) was computed by pooling the results of all 10-folds of the cross-validation72. The choice of parameters and ensemble modelling were chosen to control for potential overfitting. Moreover, we selected only three covariates for each strain distribution model to further avoid overfitting. The shape of the smoothed splines was then inspected and all had reasonable smoothing (Supplementary Fig. 3). To further improve the robustness of the approach, we performed ensemble modelling by averaging the predictions of GAMs built on each cross-validation fold using an elastic net linear regressor. Differential alpha values of 0, 0.5, and 1 were tested and the best one was retained based on the error reported by the cv.glmnet function of the glmnet R package (v4.1_7)73. Using 5-fold cross-validation, the function automatically computes the best lambda value. We avoided autocorrelation at the GFS level (as two reaches were sampled from each GFS) by randomly sampling one of the two reaches before creating each of the sub-models. By fitting nine separate models to predict each cross-validation fold, we allowed a fraction to represent local variations in the data, and stacking multiple models allowed for all samples to be included in the final model (since only one sample per GFS is included in each model). The final predictions were the mean of the predictions made by all ten ensemble models trained on 9-folds.

Response curves were obtained using predictions for all sampling points. The consistency of the stream parameter models created for each SSP scenario was assessed by comparing the selected features and their corresponding response curves (Supplementary Fig. 3). We measured model performance as cross-validation R2 comparing predicted and observed data for the held-out folds (Supplementary Table 2). Statistics comparing present and future stream parameters were obtained in the same way as statistics for the climatic dataset. Correlation between predicted stream parameters was assessed using Spearman correlations across all sites.

Models of strain and KO abundances
Models of strain abundance (normalised using the recA gene coverage) and KO abundance were built using climatic, glaciological, and mineralogical data, as well as stream parameter forecasts as covariates, chosen by feature selection (one set of features for each strain). To reduce overfitting and collinearity, we reduced the dimensionality of bioclimatic variables using principal component analysis (PCA) implemented in the prcomp R function, and retained the first six dimensions, as they represented more than 95% of the variation (Supplementary Fig. 8). Strain and KO abundances, as well as covariates were log-transformed and z-scaled as described above. Models were created with the same approach as for the stream parameters models, using GAMs stacked with an elastic net, the only difference being that spatial splines were not used; instead latitude, absolute latitude, and longitude, along with elevation and the slope of the stream, were added as potential covariates in the feature selection process. For each model, four variables were selected with the same feature selection procedure as for the stream parameters models. Model performance (predictive power) was assessed using the R2 pooled over all ten of the cross-validation folds34.

Statistical analyses
Phylogenetic signal in log2 fold-change in future abundance and in R2 prediction was measured by using the “lambda” method of the phytools R package (v1.5_1)74. We report median change and interquartile ranges for the site-specific changes for all parameters. To test for a significant shift, we conducted Wilcoxon tests of the future projections minus the present observations. Evenness was assessed by comparing present and future median predictions of abundance for all strains, and by showing that the relationship had a slope less than one (Supplementary Fig. 9). The largest monophyletic sub-trees composed only of representatives predicted to decrease under the SSP3 scenario were identified using R packages phytools (v.1.5_1)74 and ape (v5.7_1)75. The phylogenetic tree was plotted using R packages ggtree (v3.8.0)76 and ggtreeExtra (v1.10.0)77.

We used Spearman correlations between the number of shared covariates and phylogenetic distance of pairs of MAGs to assess phylogenetic signal in covariates (i.e. whether phylogenetically similar MAGs share more covariates than phylogenetically distant MAGs). Similarly, to assess taxonomic similarities, we split MAGs based on taxonomy and used the Kruskal–Wallis test to evaluate whether covariates have similar importance in all classes. This was computed using the kruskal.test R function, and plots were created to show the distribution of median relative ranks across taxa with the ggridges R package (v0.5.4)49.

Random forest classifiers were used to identify KOs associated with strains that decrease in abundance. To this end, we accounted for phylogenetic structure by splitting the data into ten phylogenetic clusters, and training models on nine of them while assessing the importance on the one left out. Random forest classifiers were created using the ranger R package (v0.15.1)78, using random grid search (n=50) to find hyperparameter settings (Supplementary Table 8). Feature importance (i.e. the importance of KOs) was estimated using the method developed by ref. 79. We considered KOs with a p-value < 0.05 and in the upper 95th percentile of the importance values in at least one of the phylogenetic clusters as significant. We considered “top” KOs that were significant in at least eight out of the ten phylogenetic clusters (n = 21), and for these, descriptions were gathered on the KEGG website (https://www.genome.jp/kegg/) (Supplementary Table 6).

Enrichment analysis was carried out at the level of KEGG categories to identify categories overrepresented in the set of significant KOs. This was done with the fisher.test in R and p-values were corrected using the Bonferroni method. We only considered positively significant categories with the thresholds: p-value < 0.05 and odds ratio (OR) > 1. To compare the number of KOs among genomes of strains predicted to decrease in abundance with “others”, we used GAMs considering the completeness, contamination, and N50 of these genomes and their interactions with a tensor (k=3, bs=’cs’). GAMs (gaussian family model) were created for the counts of KOs in the MAGs. Weights were added by multiplying the inverse of the number of MAGs for each phylogenetic cluster (such that all phylogenetic clusters obtained equal weights), the completeness of the genomes (to give less weight to KO absences owing to incomplete genomes), and the mean relative abundance. These models were fit using the bam function in the mgcv R package (v1.9_0)67. Using the same approach, we tested for differences between the two groups of bacterial genomes (i.e. “decreasing in abundance” vs “others”) in the total number of KOs, genome length and KO redundancy (unique KOs/total KOs). We further regressed the number of KOs on the genomes against genome length (accounting for completeness, contamination and N50 with a tensor; function ‘te’ with parameters k=3 and bs=’cs’). For these models, we tested the difference between the strains predicted to decrease in abundance and “others” by fitting fixed effects (and an interaction for the KO number~genomelength model) with ANOVAs as implemented in the stats R package included in the base R package (v4.3.1)46 and reported the estimated means and standard errors.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript is offering the first prediction of the change in microbial community structure of microbial communities in\nglacier fed streams under different climate change scenarios. A very timely and most interesting concept.\nI consider this manuscript groundbreaking since it is the very first approach to deal with the consequences of climate change\nand the shift of microbiome structure of glacier fed streams. Our vanishing glaciers are not only changing landscapes,\nhydrological water cycles but of course also microbial community structures in terms of composition and function.\nThe authors can successfully proof by the combination of 2.333 procaryotic metagenome-assembled genomes with climatic,\nglaciological and environmental data from glaciers worldwide that whole bacterial clades - and hence - also functionalities\nwill be lost. As a consequence, this also effects ecosystem functioning.\nIt has often been speculated that our glaciers are becoming greener with increasing temperature. Here we have a solid proof\nunder various climate change scenarios which I find extremely appealing and noteworthy.\nThis work is of significance for the field including related fields. I am not aware of any similar studies up to date. To start with,\nthere is an extreme effort in getting these sample sets together which have been collected in relatively short time in intense\nsampling campaigns (also to avoid any temporal bias [which is, however, not very likely to exclude at all considering\ndifferent seasons, weather conditions, situations after precipitation, etc.]). The project Vanishing Glaciers provided the\nground for this effort which has been successfully done in very remote areas. So far, most studies focus on one area or even\none glacier from which it is very tricky and sometimes very speculative to scale up. This manuscript compiles the major\nmountain ranges except Antarctic glaciers.\nTo my knowledge this is also the first attempt for a hierarchical, machine-learning based\nmodelling framework to forecast the\nconsequences. The model frameworks GloGEM and CHELSA provide robust base for their projections.\nDespite the fact that it is somehow a logical consequence of glaciers “to become greener” with climate change due to\nincreasing phototrophy, it is shown in the most elegant way and backed up by environmental data extremely well. The strain\ndistribution model for 2.333 strain-level MAGs is a suitable tool therefore and also explained very well (despite\nSupplementary Figure 1 is very difficult to read (text within text fields).\nWhat I also found very intriguing in this study is the prediction of a possible promotion of a parasitic lifestyle, e.g., by\nPaceibacteria due to a larger size of a host pool.\nSampling also very often offers a massive bias. Here, it has been done in a very consistent way by using the same\nprocedures and probes for environmental data. Also, considering the extreme remote character of the sampling sites, the\nsamples have been processed already in the field to avoid any change. For DOC: have samples been acidified in the field\nbefore shipping?\nThe methodology is very robust and backed up in a solid way. Metagenomic analysis are followed by state-of-the-art\nprotocols and also here they have been conducted consistently with the same procedures to avoid any bias.\nCodes are proofed and are reproducible. Generally, there is enough detail provided for the methodology to be reproduced. I\nwould not see any caveats which would require additional work.\nIn their manuscript, the authors claim the loss of clades under the aspect of climate change and hence, also the change of\necosystem functionalities. They proof it well with solid methodology, sampling and protocols. I do not see any flaws in the\nconcept as well as the presentation of the data. Figures and supplementary ones show the predicted change in the different\nvariables - also on metabolic pathways - very well.\nI am very much in favor with this study, the conceptualism thereof, the implementation,\ninterpretation and presentation.\nSupplementary information is sufficient and provides valuable background data.\nThis paper is a crucial piece of work showing us the consequences of climate change and the role of microbial communities\nwhich literally can change the functionalities in the future. I consider it as crucial to present it to decision makers for actions\nagainst climate warming. Chapeau to the authors.I highly recommend it being published.", "reviewer_2": "Bourquin et al. created a model to explore the impact that glacial-fed streams will experience under three different SSP\nclimate change scenarios. The model was built using a dataset of samples collected from a wide selection of glacial-fed\nstreams, ensuring it is based on robust data. The results show that shifts in glacial discharges due to glacial retreat will lead\nto an increase in primary productivity, mainly due to decreased turbidity, and a loss of microbial\nclades. This study is well-\npresented and clear, and it is of utmost relevance, making it suitable for publication in Nature Communications. However, I\nbelieve it could be improved by emphasizing why certain choices were made, adding some missing methodologies, and\npotentially expanding on some topics. Please find my detailed comments below:\n1. L29-32: You mention both an increase in biodiversity and a loss of entire bacterial clades under various climate change\nscenarios. These two aspects could be seen as contrasting. Could you add a brief discussion of this in your work?\n2. L85-87: Could you justify why you chose to focus on SSP1, SSP5, and particularly SSP3?\n3. L98: Is it phosphorus or phosphate? Check Supplementary Fig. 2.\n4. L113-115: While it is true that glacial flour and sediment will decrease with a decrease in glacial meltwater discharge,\nnutrient input from glaciers to the outflow systems might also decrease (for example, it looks like DIN and phosphate will\ndecrease under SSP3 and SSP5; Supplementary Fig. 2). How will this impact the “GFS greening”? Also, from\nSupplementary Fig. 2, it appears that the main change in terms of median values will be water temperature. Therefore, it\nwould be interesting to include some discussion on how this could influence algal activity.\n5. L135: Define what GAM stands for.\n6. L207: As the majority of the strains are expected to increase (L170), why did you only delve into the genomic and\nfunctional traits associated with decreased abundance? It might also be interesting to focus on the genomic traits and\nfunctional potential of strains that are expected to increase in abundance.\n7. L239-241: Wasn’t it possible to check this with your model? Maybe by looking at genes involved specifically in\nchemolithoautotrophy?\n8. L252: Should this be \"predicted to decrease\"?\n9. L293: What was the distance among the sampling patches?\n10. L322: There is a missing space.\n11. L350: I would clarify why samples from the Alaska Range were not sequenced.\n12. L375: Which database did you use when running eggNOG-mapper? You worked with the KEGG classification, not the\neggNOG classification. Is it correct?\n13. L382-383: It would be helpful to include the MAG taxonomic information (with also specifics such as completeness and\ncontamination) in the Supplementary Information (SI).\n14. I couldn’t find where you reported how chlorophyll-a and bacterial abundance were measured/obtained. Could you\nplease include this information in the Methods section?\n15. Some R packages are missing references (e.g., ggtree, ggridges).\n16. L442: Correct “as well as and covariates” to “as well as covariates.”\n17. Fig 2: Do the “Group” colors correspond to the branch colors in the phylogenetic tree? Why are “not significant” and\n“increase” strains clustered and represented together as just one group? I would also suggest flagging the clades predicted\nto increase with a different color, as this is also an important result.\n18. Supplementary Fig. 1: Please specify what the different colors in panel B represent.", "reviewer_3": "This is one of the first manuscripts that I've encountered that attempted to examine the likely changes to environmental\nmicrobial community structure and function under different climate scenarios. They specifically focused on glacier fed\nstreams (GFS)--an important environment as glacier meltwater is a significant contributor to freshwater resources globally\nand is rapidly changing as a result of melt due to climate change. The work is highly significant, as generalized climate\nmodels usually predict the effects of a warming atmosphere on physical parameters (e.g. future climate, ocean circulation,\netc.), chemical parameters (e.g. weathering rates and chemical products, carbon and nutrient cycling, etc.), and biological\nparameters (e.g. plant growth and coverage, species shifts, etc.); however, very little has been done to predict the impact of\nclimate change on specific microbial communities. Since these communities underlie many of the important ecosystem\nprocesses and provide significant feedbacks to the climate system, it is crucial to understand how they will respond to\nenvironmental change.\nThe general approach used in this manuscript was to use a space-for-time approach where one sample was collected close\nto the glacier terminus (\"early\") and another downstream (\"late\"). They then examined the differences between the \"early\"\nand \"late\" samples in terms of biogeochemistry, biomass, and alpha diversity metrics. Using Generalized Additive\nModels (GAMs), they determined the correlation of these parameters with climate change, and predicted how they would change in\nthe future under three climate models (SSPs). These parameters were correlated with a large set of strain-resolved\nmetagenome assembled genomes (MAGs) to examine strain responses to the predicted changes as well as examining the\nresponse of functions (based on predicted changes in KEGG orthologs and KEGG categories). They argue that cross\nvalidation from several different predictive approaches (combination of various GAMs and random forest models) makes\ntheir predictive approach robust.\nAltogether, they predict 1) a shift toward more generalist and away from more specialist strains, 2) a shift toward light\nacquisition pathways, which is consistent with a future \"greening\" of GFS). This could lead to the loss of entire clades and a\nmajor shift in carbon and nutrient cycling in these systems.\nOverall, I found this paper and their approach exciting and their results fascinating. I have a few concerns, some significant\nand others less so.\n1. Their space-for-time approach is not well explained. The approach taken might well be reasonable, but it is unclear how\nthe conversion of space for time was accomplished, especially considering that a) there are only two points at each location\nfor space, b) glaciers retreat at dramatically different rates, c) the sediments that they sampled (sandy sediments) are quite\nmobile in GFS, and may not be the same age as the location (and this would be affected by the elevation changes in the\nsystem) and d) the linearity of the space-for-time assumptions is not established. Some further description and justification of\nthe space-for-time approach and its validity would be welcomed.\n2. If we accept that their models are correct, they make predictions regarding shifts in mean nearest taxon distance, strain\ndistribution, genomic traits, and functional potential. In the cases where these features can be directly correlated to the\nenvironmental parameters, these predictions are reasonable. However, in some cases they go well beyond that direct\ncorrelation. For example, they predict that Gamma and Alpha Proteobacteria will increase in abundance, which is\nreasonable based on the correlation of strains from these clades with environmental parameters. However, they then go on\nto argue that \"Gamma- and Alphaproteobacteria in GFS are also associated with thiamine (vitamin B1) and vitamin B12\nbiosynthesis, and putatively can alleviate the demands of auxotrophic algae\". This argument is based on a broad correlation\nof these phylum-level clades with those features, and do not specifically relate to the features of the strains characterized in\nthis study. This happens several times in the main text (e.g. lines 181-185, 220-222, 254-255). Only proximal and direct\npredictions based upon specific changes in their models should be included. Don't go too far beyond the data.\n3. While predictions of environmental parameters and biomass based on the models, it is not at all clear to me how alpha\ndiversity is being predicted here based on the environmental data (lines 120-124). The description in the methods is simply\n\"Biomass and diversity metrics for the bacterial communities were modelled analogously.\" (lines 389-390). Could you\nexplain this more clearly?\n4. For the GAMs, which samples were used to train the models and which to test them? This is not clear. Please clearly\nstate.\n5. There was no metagenomics performed on the Alaska Range (line 350). Why? Please clarify.\n6. I didn't quite follow the co-assembly approach for constructing their MAGs. They state \"For each individual assembly, we\nmapped the reads of the 5closest samples (Euclidean distances of gps coordinates)\" (lines 362-363). Are they co-\nassembling MAGs from the same location (i.e. both reaches along the stream) as well as the four other closest locations (by\nEuclidean distance), or is this something else? Why did they approach your co-assembly this way--seems rather random.\n7. They don't describe variability very much. While I understand that they are attempting to make predictions about GFS\nresponse globally to climate change, I imagine that there must be some fairly signficant variation depending upon location\nas well as other parameters. This will likely require one or several other manuscripts to describe in detail; however, some\ncharacterization of how representative their results are and how likely any one particular GFS is likely to see these same\nchanges would help strengthen their argument here.\n8. Similar to point 7, the sole or at least primary driver described here for GFS benthic communities is glacial melt input.\nWhile this is true near the glacier, as one gets further from the source, the influence of other components (e.g. plant organic\nmatter, sediment and soil erosion, meteorological inputs, anthropogenic inputs) increases. So while their predictions may\nhold for GFS close to the glaciers, it is unclear to me how well this will translate to changes further downstream. This also\nadds a bit to my concerns regarding their space-for-time substitution approach. Some statements (caveats or clarifications)\nregarding this would be helpful.\nOnce again, I am excited and broadly convinced regarding their approach and their findings. Their findings are consistent\nwith other predicted changes resulting from glacial melt, but makes the shifts quantitative and palpable. It is generally well\nwritten and easy to follow. I applaud the authors for making a very complex paper very accessible. I think this paper will be\nused as a guide for many future studies. Kudos.", "reviewer_author_response": "Response to the reviewers for the manuscript “Predicting climate-change\nimpacts on the global glacier-fed stream microbiome”\nReviewer #1\nThis manuscript is offering the first prediction of the change in microbial community structure\nof microbial communities in glacier fed streams under different climate change scenarios. A\nvery timely and most interesting concept.\nI consider this manuscript groundbreaking since it is the very first approach to deal with the\nconsequences of climate change and the shift of microbiome structure of glacier fed\nstreams. Our vanishing glaciers are not only changing landscapes, hydrological water cycles\nbut of course also microbial community structures in terms of composition and function.\nAUTHORS: We are delighted to receive such positive feedback to our work. We would like\nto express our sincere gratitude to reviewer #1 for these comments.\nThe authors can successfully proof by the combination of 2.333 procaryotic metagenome-\nassembled genomes with climatic, glaciological and environmental data from glaciers\nworldwide that whole bacterial clades - and hence - also functionalities will be lost. As a\nconsequence, this also effects ecosystem functioning.\nIt has often been speculated that our glaciers are becoming greener with increasing\ntemperature. Here we have a solid proof under various climate change scenarios which I find\nextremely appealing and noteworthy.\nThis work is of significance for the field including related fields. I am not aware of any similar\nstudies up to date. To start with, there is an extreme effort in getting these sample sets\ntogether which have been collected in relatively short time in intense sampling campaigns\n(also to avoid any temporal bias [which is, however, not very likely to exclude at all\nconsidering different seasons, weather conditions, situations after precipitation, etc.]). The\nproject Vanishing Glaciers provided the ground for this effort which has been successfully\ndone in very remote areas. So far, most studies focus on one area or even one glacier from\nwhich it is very tricky and sometimes very speculative to scale up. This manuscript compiles\nthe major mountain ranges except Antarctic glaciers.\nAUTHORS: We are grateful to this reviewer for recognizing the efforts that went into this\nstudy — ranging from the field work to sequencing and modeling.\nTo my knowledge this is also the first attempt for a hierarchical, machine-learning based\nmodelling framework to forecast the consequences. The model frameworks GloGEM and\nCHELSA provide robust base for their projections.\nDespite the fact that it is somehow a logical consequence of glaciers “to become greener”\nwith climate change due to increasing phototrophy, it is shown in the most elegant way and\nbacked up by environmental data extremely well. The strain distribution model for 2.333\nstrain-level MAGs is a suitable tool therefore and also explained very well (despite\nSupplementary Figure 1 is very difficult to read (text within text fields).\nAUTHORS: We are very happy to read this positive appraisal of our modeling framework.\nWe do agree that a hierarchical structure, building on robust projections such as those\nprovided by GloGEM and CHELSA, makes sense for predictive modeling of microbial\npopulations. Following this comment, we have revised Supplementary Figure 1 to increase\nreadability and clarity.\nparasitic lifestyle, e.g., by Paceibacteria due to a larger size of a host pool.\nSampling also\nvery often offers a massive bias. Here, it has been done in a very consistent way by using\nthe same procedures and probes for environmental data. Also, considering the extreme\nremote character of the sampling sites, the samples have been processed already in the\nfield to avoid any change. For DOC: have samples been acidified in the field before\nshipping?\nAUTHORS: We agree with this reviewer on the necessity to properly handle the various\nsamples — starting right in the field. Therefore, we developed SOPs for the Vanishing\nGlacier Project that were strictly and consistently applied across all study sites as variation in\nsampling procedures, sequencing and data processing can inflate variance of microbial\ncommunities, making predictions very hard.\nRegarding DOC samples - they were filtered (ashed GF/F filters) and stored in pre-\ncombusted glass vials at 4°C without acidification. Since shipping of acids is very\nproblematic, we would have had to rely on local supplies of acid of the highest grade, which\nunfortunately is often not (or not easily) available in the countries that we visited. Therefore,\nto avoid potential contamination, we preferred not to acidify our DOC samples in the field.\nHowever, we have been running assays to determine the biologically available fraction of\nDOC (BDOC) using Alpine GFS samples. For this, filtered DOC samples are amended with\na native bacterial community, and DOC concentration is monitored over 28 days\n(incubations in the dark and at room temperature). Typically, we find very low BDOC,\nranging between 0 and 20% of initial DOC concentrations. We therefore consider the error\nintroduced due to storage of our DOC samples to have minor influence.\nThe methodology is very robust and backed up in a solid way. Metagenomic analysis are\nfollowed by state-of-the-art protocols and also here they have been conducted consistently\nwith the same procedures to avoid any bias.\nAUTHORS: Thank you for this comment, which we fully agree with.\nCodes are proofed and are reproducible. Generally, there is enough detail provided for the\nmethodology to be reproduced. I would not see any caveats which would require additional\nwork.\nAUTHORS: This is indeed an important aspect of our study and we are grateful to the\nreviewer for highlighting it. We hope that the availability of our code helps other researchers\nadopt a similar, hierarchical modeling framework.\nIn their manuscript, the authors claim the loss of clades under the aspect of climate change\nand hence, also the change of ecosystem functionalities. They proof it well with solid\nmethodology, sampling and protocols. I do not see any flaws in the concept as well as the\npresentation of the data. Figures and supplementary ones show the predicted change in the\ndifferent variables - also on metabolic pathways - very well.\nAUTHORS: We are grateful for this positive tone here as well.\nI am very much in favor with this study, the conceptualism thereof, the implementation,\ninterpretation and presentation. Supplementary information is sufficient and provides\nvaluable background data.\nthe role of microbial communities which literally can change the functionalities in the future. I\nconsider it as crucial to present it to decision makers for actions against climate warming.\nChapeau to the authors. I highly recommend it being published.\nAUTHORS: Overall, we would like to reiterate our gratitude for this reviewer for her/his\npositive feedback and the recognition of our work.\nReviewer #2\nBourquin et al. created a model to explore the impact that glacial-fed streams will experience\nunder three different SSP climate change scenarios. The model was built using a dataset of\nsamples collected from a wide selection of glacial-fed streams, ensuring it is based on robust\ndata. The results show that shifts in glacial discharges due to glacial retreat will lead to an\nincrease in primary productivity, mainly due to decreased turbidity, and a loss of microbial\nclades. This study is well-presented and clear, and it is of utmost relevance, making it\nsuitable for publication in Nature Communications. However, I believe it could be improved\nby emphasizing why certain choices were made, adding some missing methodologies, and\npotentially expanding on some topics. Please find my detailed comments below:\nAUTHORS: We are pleased to read that reviewer # 2 finds our work to be of great\nrelevance. We thank the reviewer for the constructive feedback, which we address below.\n1. L29-32: You mention both an increase in biodiversity and a loss of entire bacterial clades\nunder various climate change scenarios. These two aspects could be seen as contrasting.\nCould you add a brief discussion of this in your work?\nAUTHORS: Thank you for pointing this out — this is an important comment, and we\napologize for not having been clear enough on this. Indeed, bacterial biodiversity (as\nShannon H’) is predicted to increase by ~ 6% in future GFSs. We speculate that this is\nexplicable by increased niche availability linked to increased resource availability (due to\n“greening”). By the way, an accompanying paper, accepted in Nature Microbiology (Michoud\net al.), suggests indeed that biofilms in GFSs with elevated algal biomass are more complex,\nhence potentially offering more diverse niches. Furthermore, present-day GFS microbial\ndiversity is dominated by microbes that are well adapted to the selective constraints of the\nGFS environment. For instance, we previously linked environmental constraints to enhanced\nmicrodiversity in GFS bacterial communities (doi.org/10.1038/s41396-021-01106-6).\nMicrodiversity (i.e. the presence of many closely-related taxa) can be assessed using\nphylogenetic measures, such as nearest taxon distance. Low nearest-taxon-distances reflect\nphylogenetic clustering - which is high in communities under selection. In future GFSs,\nincreased resource availability may weaken these selective constraints and allow for more\nbiodiversity, however, at the expense of these well-adapted and microdiverse clades. Our\nmodel predictions are in line with these expectations - showing an increase in nearest taxon\ndistances (reflecting the loss of microdiversity), the putative loss of entire clades, and at the\nsame time higher overall diversity. We revised the text to better explain the apparent\ndiscrepancy between a loss of clades and a concomitant increase in biodiversity.\n2. L85-87: Could you justify why you chose to focus on SSP1, SSP5, and particularly SSP3?\nAUTHORS: These simply reflect extreme positive (SSP1, also called “Sustainability - Taking\nthe Green Road”), extreme negative (SSP5, also called “Fossil-fueled Development - Taking\n3\nthe Highway”) and the intermediate (SSP3, also called “Regional Rivalry - A Rocky Road”)\nclimate change scenarios. These scenarios are often used in predictive studies in ecology.\n3. L98: Is it phosphorus or phosphate? Check Supplementary Fig. 2.\nAUTHORS: Thanks for noticing this - phosphorus concentration in the streamwater was\nmeasured as soluble reactive phosphorus (SRP). We corrected this in Supplementary Fig. 2.\n4. L113-115: While it is true that glacial flour and sediment will decrease with a decrease in\nglacial meltwater discharge, nutrient input from glaciers to the outflow systems might also\ndecrease (for example, it looks like DIN and phosphate will decrease under SSP3 and\nSSP5; Supplementary Fig. 2). How will this impact the “GFS greening”? Also, from\nSupplementary Fig. 2, it appears that the main change in terms of median values will be\nwater temperature. Therefore, it would be interesting to include some discussion on how this\ncould influence algal activity.\nAUTHORS: This is an interesting comment for which are very grateful. To be honest, and\nthis is the beauty of the science that we all do, we do not know. We discussed this issue to\nsome extent already in a previous paper (doi.org/10.1038/s41561-024-01393-6). Our\nmodeling framework forecasts both a decrease in turbidity and concomitantly, predicts\nincreases in DIN and SRP. We acknowledge that GFS nutrient concentration depends on a\nmultitude of factors, including underlying geology, atmospheric deposition, and catchment\nvegetation and related soil development. For instance, climate-induced glacier shrinkage\nreduces the erosive capacity of glaciers, thereby potentially impacting phosphorus supply\nfrom bedrock. Smaller glaciers may also receive reduced atmospheric deposition\n(particularly affecting N in GFS). However, our previous study (doi.org/10.1038/s41561-024-\n01393-6; Fig. 1b) across a wide range of GFSs showed that benthic algal biomass is\nstronger linked to turbidity (i.e., hence light availability) than to the availability of inorganic\nnutrients. Hence, our reflections on a future “greening” of GFSs and streamwater turbidity\nseem supported.\n5. L135: Define what GAM stands for.\nAUTHORS: Thank you for pointing this out. GAM stands for Generalized Additive Model. We\ncorrected this in the revised manuscript.\n6. L207: As the majority of the strains are expected to increase (L170), why did you only\ndelve into the genomic and functional traits associated with decreased abundance? It might\nalso be interesting to focus on the genomic traits and functional potential of strains that are\nexpected to increase in abundance.\nAUTHORS: This is a valid point for which we are grateful. Our choice was motivated by the\nfact that we have genomic information available for microbes abundant in present-day GFSs,\nincluding those that are predicted to decrease in abundance. Focusing on these strains\ntherefore provides robust interpretation.\n7. L239-241: Wasn’t it possible to check this with your model? Maybe by looking at genes\ninvolved specifically in chemolithoautotrophy?\nAUTHORS: No, it is not possible to obtain this information directly from our models. The\nmodels predict shifts in the abundance of MAGs and from these we then inferred changes in\nthe functional potential of the microbiome by looking at the genes that the respective MAGs\nencode. Directly resolving individual genes that may be more or less prevalent in the future,\nwould be valid and interesting but would (potentially) require a different modeling framework\n(e.g considering functional redundancy, orthologs, etc).\n8. L252: Should this be \"predicted to decrease\"?\nAUTHORS: Yes. We corrected this mistake.\n9. L293: What was the distance among the sampling patches?\nAUTHORS: Sediment patches were sampled within less than 10m apart. We added this\ninformation to the revised manuscript.\n10. L322: There is a missing space.\nAUTHORS: Thank you. We corrected this mistake.\n11. L350: I would clarify why samples from the Alaska Range were not sequenced.\nAUTHORS: Alaska was our last expedition and samples simply came in too late for\nsequencing. However, we opted to maintain the Alaska GFSs for modeling the environment.\n12. L375: Which database did you use when running eggNOG-mapper? You worked with\nthe KEGG classification, not the eggNOG classification. Is it correct?\nAUTHORS: We used eggNOG-mapper (v2.1.9) against the eggNOG v5.0 database, which\nalso contains KEGG orthologues. We added the database to the revised manuscript.\n13. L382-383: It would be helpful to include the MAG taxonomic information (with also\nspecifics such as completeness and contamination) in the Supplementary Information (SI).\nAUTHORS: We refer the reviewer to the accompanying paper in Nature Microbiology for a\ndetailed breakdown of MAG-related information.\n14. I couldn’t find where you reported how chlorophyll-a and bacterial abundance were\nmeasured/obtained. Could you please include this information in the Methods section?\nAUTHORS: Thank you for pointing this out. We apologize for this lapsus and have now\nadded the respective information in the Methods section.\n15. Some R packages are missing references (e.g., ggtree, ggridges).\nAUTHORS: Thank you for spotting this, we have included the missing references.\n16. L442: Correct “as well as and covariates” to “as well as covariates.”\nAUTHORS: Corrected.\n17. Fig 2: Do the “Group” colors correspond to the branch colors in the phylogenetic tree?\nWhy are “not significant” and “increase” strains clustered and represented together as just\none group? I would also suggest flagging the clades predicted to increase with a different\ncolor, as this is also an important result.\nAUTHORS: Thank you for pointing this out. Given that we focused on MAGs predicted to\ndecrease in abundance, we initially represented only the decreasing clades on this\nphylogenetic tree. We have now added clades predicted to increase in the tree.\n18. Supplementary Fig. 1: Please specify what the different colors in panel B represent.\ncomment by Reviewer # 1) and clarified the meaning of the different colors (i.e. dataset\nsources) in the legend.\nReviewer #3\nThis is one of the first manuscripts that I've encountered that attempted to examine the likely\nchanges to environmental microbial community structure and function under different climate\nscenarios. They specifically focused on glacier fed streams (GFS)--an important\nenvironment as glacier meltwater is a significant contributor to freshwater resources globally\nand is rapidly changing as a result of melt due to climate change. The work is highly\nsignificant, as generalized climate models usually predict the effects of a warming\natmosphere on physical parameters (e.g. future climate, ocean circulation, etc.), chemical\nparameters (e.g. weathering rates and chemical products, carbon and nutrient cycling, etc.),\nand biological parameters (e.g. plant growth and coverage, species shifts, etc.); however,\nvery little has been done to predict the impact of climate change on specific microbial\ncommunities. Since these communities underlie many of the important ecosystem processes\nand provide significant feedbacks to the climate system, it is crucial to understand how they\nwill respond to environmental change.\nAUTHORS: We would like to express our gratitude also to this reviewer for recognizing the\nsignificance of our work. Quite obviously, we strongly agree with the reviewer regarding the\nimportance of these unique microbial communities threatened by climate change.\nThe general approach used in this manuscript was to use a space-for-time approach where\none sample was collected close to the glacier terminus (\"early\") and another downstream\n(\"late\"). They then examined the differences between the \"early\" and \"late\" samples in terms\nof biogeochemistry, biomass, and alpha diversity metrics. Using Generalized Additive\nModels (GAMs), they determined the correlation of these parameters with climate change,\nand predicted how they would change in the future under three climate models (SSPs).\nThese parameters were correlated with a large set of strain-resolved metagenome\nassembled genomes (MAGs) to examine strain responses to the predicted changes as well\nas examining the response of functions (based on predicted changes in KEGG orthologs\nand KEGG categories). They argue that cross validation from several different predictive\napproaches (combination of various GAMs and random forest models) makes their\npredictive approach robust.\nAUTHORS: This summary of our modeling framework is correct - and we appreciate the\nmore detailed comments below which, in our opinion, further strengthened the manuscript.\nAltogether, they predict 1) a shift toward more generalist and away from more specialist\nstrains, 2) a shift toward light acquisition pathways, which is consistent with a future\n\"greening\" of GFS). This could lead to the loss of entire clades and a major shift in carbon\nand nutrient cycling in these systems. Overall, I found this paper and their approach exciting\nand their results fascinating. I have a few concerns, some significant and others less so.\nAUTHORS: Thank you again for the positive and constructive feedback.\n1. Their space-for-time approach is not well explained. The approach taken might well be\nreasonable, but it is unclear how the conversion of space for time was accomplished,\n6\nespecially considering that a) there are only two points at each location for space, b) glaciers\nretreat at dramatically different rates, c) the sediments that they sampled (sandy sediments)\nare quite mobile in GFS, and may not be the same age as the location (and this would be\naffected by the elevation changes in the system) and d) the linearity of the space-for-time\nassumptions is not established. Some further description and justification of the space-for-\ntime approach and its validity would be welcomed.\nAUTHORS: We thank the reviewer for these thoughtful comments. We extended the\nmanuscript to better address these points. Space-for-time substitution is a powerful\napproach to study ecological processes over longer time-scales and has been repeatedly\nused in ecological studies (e.g., doi.org/10.1038/s41559-017-0426-x;). However, we agree\nthat certain criteria require careful consideration for this approach. Importantly, work on\nvegetation succession in proglacial forefields (doi/abs/10.2307/2937039) has shown that\nspace-for-time substitution is particularly useful when deterministic processes drive the\nassembly of rather low-diversity communities where locations along the chronosequence\nfollow similar ecological trajectories doi.org/10.1111/j.1365-2745.2010.01664.x).\nIn this context, our previous work (doi.org/10.1038/s41396-021-01106-6) showed that\nbenthic bacterial communities in GFS are indeed characterized by a rather low diversity and\nshaped by deterministic assembly (i.e. homogeneous selection). Our selection of sampling\npoints was motivated by the necessity to maximize the space-for-time substitution gradient -\nbut at the same time constrained by being part of the same ecological trajectories. We\ntherefore sampled sediments from an upstream reach as close as possible to the glacier\nsnout (typically a few meters from the snout; limited only by field safety aspects) and a reach\nclose to the terminal moraine demarking the last glacial maximum (i.e., the “Little Ice Age”;\nLIA) approximately 150 to 170 years ago. Below LIA, soils are typically better developed and\nstream channels more stable. This induces a rather non-linear change in environmental\nconditions. However, it is reasonable to assume at least a quasi-linearity upstream from LIA.\nWe acknowledge that this is complex, depending on glacier length recession dynamics (also\nincluding periods of positive mass balances). It is unfortunate that this cannot simply not be\nresolved for each glacier individually, at least not within the scope of this study. We also\nwant to draw the attention to the fact that we were simply limited logistically to two reaches\nper GFS only. On average, sampling one single reach took up to four hours, hence the\ntriplicate time per GFS would simply not have been feasible. However, we used penalized\nsmooth splines (i.e. non-linear relationships smoothed with a penalty on the number of\nknots) to model relationships between covariates.\nRegarding the sandy sediment fraction - it is correct that it is constantly transported\ndownstreams in GFS therefore potentially reducing the differences between up- and\ndownstream benthic microbial communities. Bacterial dispersing along their sedimentary\nsubstrate may indeed add random noise to these communities - nevertheless we find that\ndeterministic assembly processes predominantly shape benthic bacterial communities in\nGFS. Moreover, given the importance of algal biomass highlighted in previous work\n(doi.org/10.1038/s41561-024-01393-6) and our models, we argue that the time-scales of\nsediment reworking (in downstream reaches) and growth of benthic algae appear sufficiently\nlong to sustain microbial niches. We have implemented these considerations in the revised\nmanuscript.\n7\n2. If we accept that their models are correct, they make predictions regarding shifts in mean\nnearest taxon distance, strain distribution, genomic traits, and functional potential. In the\ncases where these features can be directly correlated to the environmental parameters,\nthese predictions are reasonable. However, in some cases they go well beyond that direct\ncorrelation. For example, they predict that Gamma and Alpha Proteobacteria will increase in\nabundance, which is reasonable based on the correlation of strains from these clades with\nenvironmental parameters. However, they then go on to argue that \"Gamma- and\nAlphaproteobacteria in GFS are also associated with thiamine (vitamin B1) and vitamin B12\nbiosynthesis, and putatively can alleviate the demands of auxotrophic algae\". This argument\nis based on a broad correlation of these phylum-level clades with those features, and do not\nspecifically relate to the features of the strains characterized in this study. This happens\nseveral times in the main text (e.g. lines 181-185, 220-222, 254-255). Only proximal and\ndirect predictions based upon specific changes in their models should be included. Don't go\ntoo far beyond the data.\nAUTHORS: We are grateful for this critical comment, based on which we have removed the\nspeculation regarding interactions between bacteria and algae. However, we also believe\nthat hypotheses that arise from our models in combination with independent lines of\nevidence are interesting for other readers - but agree that they needed to be more clearly\nphrased as such. We therefore rephrased these instances to more clearly highlight them as\nhypotheses that, for instance, may be addressed in future work.\n3. While predictions of environmental parameters and biomass based on the models, it is not\nat all clear to me how alpha diversity is being predicted here based on the environmental\ndata (lines 120-124). The description in the methods is simply \"Biomass and diversity\nmetrics for the bacterial communities were modelled analogously.\" (lines 389-390). Could\nyou explain this more clearly?\nAUTHORS: Thank you for pointing out this lack of clarity. What we meant is that these\nmodels were created the same way as the environmental parameters ones. We revised this\nsection to be clearer.\n4. For the GAMs, which samples were used to train the models and which to test them? This\nis not clear. Please clearly state.\nAUTHORS: This work is based on cross-validation. Models were created and validated\nusing 10-fold cross-validation (each time training on 90% of the glacier-fed streams, and\ntesting on the remaining ones). We rephrased the respective section to better clarify this\nimportant point.\n5. There was no metagenomics performed on the Alaska Range (line 350). Why? Please\nclarify.\nAUTHORS: This point was also raised by reviewer #2. Alaska was our last expedition and\nsamples simply came in too late for sequencing. However, we opted to maintain the Alaska\nGFSs for modeling the environment.\n6. I didn't quite follow the co-assembly approach for constructing their MAGs. They state\n\"For each individual assembly, we mapped the reads of the 5 closest samples (Euclidean\ndistances of gps coordinates)\" (lines 362-363). Are they co-assembling MAGs from the same\nlocation (i.e. both reaches along the stream) as well as the four other closest locations (by\nEuclidean distance), or is this something else? Why did they approach your co-assembly this\nway--seems rather random.\nAUTHORS: We apologize for this unclarity. In fact, we assembled reads into contigs from\neach GFS reach individually (i.e. technically spoken they were not co-assembled). Then, to\nbin contigs into MAGs, we mapped reads from five samples to the contigs. This is called\nmulti-coverage binning, which is now commonly done and increases the completeness and\ndecreases contamination of MAGs (doi.org/10.1038/s41592-023-01934-8). Unfortunately,\nmulti-coverage binning is computationally expensive and mapping reads from all\nmetagenomes was infeasible. Hence, we mapped reads of five geographically near sites,\nincluding up- and downstream reaches of the same GFS (where available). We revised this\nsection to be clear about this important step.\n7. They don't describe variability very much. While I understand that they are attempting to\nmake predictions about GFS response globally to climate change, I imagine that there must\nbe some fairly signficant variation depending upon location as well as other parameters. This\nwill likely require one or several other manuscripts to describe in detail; however, some\ncharacterization of how representative their results are and how likely any one particular\nGFS is likely to see these same changes would help strengthen their argument here.\nAUTHORS: We are grateful for this very valid point. Please, note that we report median\nvalues of change across all GFSs, and where applicable interquartile ranges as well. Hence\nthis provides a good first appreciation of variability related to the global data set. As correctly\nrecognized by this reviewer, exploring forecasts further for regionality would be great but\ncertainly beyond the scope of this manuscript and not trivial either (e.g., for reasons of\nstatistical power). We have provided forecasts for three climate change scenarios. These\nforecasts do not differ qualitatively, but rather in the magnitude of individual responses.\nWhile we agree that better understanding regional (or even local) variability will be\ninteresting and relevant, this putatively requires temporal datasets (of specific locations) and\nis beyond the scope of this work.\n8. Similar to point 7, the sole or at least primary driver described here for GFS benthic\ncommunities is glacial melt input. While this is true near the glacier, as one gets further from\nthe source, the influence of other components (e.g. plant organic matter, sediment and soil\nerosion, meteorological inputs, anthropogenic inputs) increases. So while their predictions\nmay hold for GFS close to the glaciers, it is unclear to me how well this will translate to\nchanges further downstream. This also adds a bit to my concerns regarding their space-for-\ntime substitution approach. Some statements (caveats or clarifications) regarding this would\nbe helpful.\nAUTHORS: Thank you for raising this important point. Indeed, our predictions are restricted\nto the GFS reaches that are relatively close to the glaciers (within the LIA chronosequence);\nwe have made this point clearer in the revised manuscript. However, we also argue that\nbecause of the absence of major confounding variation within the upper reaches under direct\ninfluence of glaciers, GFSs represent a suitable environment for such predictions, allowing\nus to model proximal climate change-related shifts in microbial communities. We consider\nthis as a first step before downstream confounding variation (e.g. land use) can be\nconsidered in even more complex models. The reviewer is also right in the sense that our\nspace-for-time substitution approach was tailored to resolve the most recent deglaciation\nhistory - and hence allowing for a “mirrored” forecasting horizon of a few decades.\nTheir findings are consistent with other predicted changes resulting from glacial melt, but\nmakes the shifts quantitative and palpable. It is generally well written and easy to follow. I\napplaud the authors for making a very complex paper very accessible. I think this paper will be\nused as a guide for many future studies. Kudos.\nAUTHORS: We are very pleased with this positive assessment. Also we would like to\nreiterate our gratitude to this reviewer for her/his thoughtful comments. They have further\nimproved the manuscript.\n10"}, "author_response": ""}, {"round": 2, "reviewer_comments": {"reviewer_2": "The authors have thoroughly addressed my comments and suggestions, and I would like to thank them for their effort and\nattention to detail. They have made the necessary revisions to clarify key points, improve the structure, and incorporate\nadditional information where needed. In my opinion, the changes have significantly strengthened the overall quality and\ncoherence of the work.", "reviewer_3": "Great job responding to the concerns and critiques of me and the other reviewers. I am happy with the current version and\nwith the responses."}, "author_response": ""}]

Decision letter