Microbial life in glacier-fed streams (GFSs) is dominated by benthic biofilms which fulfill critical ecosystem processes. However, it remains unclear how the bacterial communities of these biofilms assemble in stream ecosystems characterized by rapid turnover of benthic habitats and high suspended sediment loads. Using16S rRNA gene amplicon sequence data collected from 54 GFSs across the Himalayas, European Alps, and Scandinavian Mountains, we found that benthic biofilms harbor bacterial communities that are distinct from the bacterial assemblages suspended in the streamwater. Our data showed a decrease in species richness in the benthic biofilms compared to the bacterial cells putatively free-living in the water. The benthic biofilms also differed from the suspended water fractions in terms of community composition. Differential abundance analyses highlighted bacterial families that were specific to the benthic biofilms and the suspended assemblages. Notably, source-sink models suggested that the benthic biofilm communities are not simply a subset of the suspended assemblages. Rather, we found evidence that deterministic processes (e.g., species sorting) shape the benthic biofilm communities. This is unexpected given the high vertical mixing of water and contained bacterial cells in GFSs and further highlights the benthic biofilm mode of life as one that is determined through niche-related processes. Our findings therefore reveal a “native” benthic biofilm community in an ecosystem that is currently threatened by climate-induced glacier shrinkage.
IMPORTANCE Benthic biofilms represent the dominant form of life in glacier-fed streams. However, it remains unclear how bacterial communities within these biofilms assemble. Our findings from glacier-fed streams from three major mountain ranges across the Himalayas, the European Alps and the Scandinavian Mountains reveal a bacterial community associated with benthic biofilms that is distinct from the assemblage in the overlying streamwater. Our analyses suggest that selection is the underlying process to this differentiation. This is unexpected given that bacterial cells that are freely living or attached to the abundant sediment particles suspended in the water continuously mix with the benthic biofilms. The latter colonize loose sediments that are subject to high turnover owing to the forces of the water flow. Our research unravels the existence of a microbiome specific to benthic biofilms in glacier-fed streams, now under major threats due to global warming.


Stream and river ecosystems are characterized by the continuous downstream flow of water. This directed force influences the spatial dynamics of biodiversity and ecological communities both at the scale of individual streams or rivers (i.e., linear) and the networks (i.e., dendritic) that they form (1). Owing to their small size, microorganisms are easily dispersed from terrestrial sources into headwaters and further downstream through fluvial networks. It is a common understanding that the dispersal (e.g., through mass effects) of microbial cells shapes microbial assemblages in the water of headwater streams while species sorting gains relevance downstream (27). Most of these studies have focused on the longitudinal (i.e., downstream) dimension of microbial community dynamics in streams and rivers. However, relatively few have attempted to unravel the vertical dimension of dispersal in streams and rivers, understood as the downwards transport of cells through the water column toward the bed sediments. Given the continuous turbulence-induced mixing of microorganisms in the water with those colonizing the benthic sediments, particularly in low-submergent (i.e., shallow water but high bed roughness) streams, it is plausible to assume homogeneity between both assemblages. However, as shown by a few studies (8, 9) benthic biofilms at the interface between the bed sediments and the overlying water may indeed select certain taxa over others.
Glacier-fed streams (GFSs) are extreme headwaters with highly unconsolidated sediments and unstable channels, and therefore exhibit frequent and pronounced turnover of their sedimentary habitats. Furthermore, high loads of fine suspended sediments in transport by the streamflow frequently mix with the benthic sediments. These hydraulic and sedimentary processes are particularly pronounced during glacier melt in summer when high discharge mobilizes sediments from various sources within the catchment, including subglacial and supraglacial environments as well as proglacial channel banks (1012). Despite this environmental harshness, GFSs support diverse life (1315). For instance, microbial biofilms that colonize various sedimentary habitats, ranging from sand to boulders, can host relatively high levels of biodiversity spanning the entire tree of life (8, 1619). Biofilms fulfill fundamental ecosystem processes in GFSs, including nutrient cycling and primary production (17, 20), and form the basis of the food chain, further supporting invertebrate communities (15). Given the high turnover of GFS sedimentary habitats, it seems plausible that dispersal of microbial cells contained in the streamwater is critical in shaping benthic biofilm communities in GFSs and “rescuing” them after disturbances. At the same time, the continuous mixing of suspended and bed sediments, and their particle-attached bacteria, might prevent benthic biofilms from establishing specific bacterial communities because of overwhelming mass effects. Metacommunity theory (21) offers a framework to assess the relative roles of dispersal and selection in shaping biofilm communities. Yet, the relative importance of these processes for community assembly in GFS benthic biofilms remains inadequately understood. Increasing our understanding about these mechanisms will be important given the rapid pace of glacier shrinkage worldwide (22) and the potential deleterious effects on GFS microbial diversity.
The aim of the present study was to explore the “vertical” dimension of microbial metacommunity dynamics across a broad range of GFSs in the European Alps, the Scandinavian Mountains, and the Himalayas. Specifically, we tested whether biofilms associated with benthic sediments in GFSs form distinct communities that assemble from the pool of bacteria suspended in the streamwater. Rather than investigating biofilms associated with boulders or rocks as previously done (8, 23), we focused on those associated with the sandy fraction of the benthic sediment. Because of their high surface-to-volume ratio, sandy sediments greatly contribute to the colonizable surface area in streams but are prone to flow-induced erosion. Using 16S rRNA gene sequencing data, we compared biodiversity and community composition patterns between benthic biofilms and their potential source assemblages in the overlying streamwater in GFSs. Next, we referred to metacommunity theory (2426) to assess the relevance of deterministic (i.e., species sorting) versus stochastic (i.e., dispersal) processes that shape the community composition of the benthic biofilms. In clear water streams and rivers, free-living cells may disproportionately form the source assemblage; given the high sediment load in GFSs, we considered both unattached (UC; putatively free-living) microorganisms and cells attached (AC) to fine suspended sediments as our source material. We anticipated that high sediment mixing and turnover in GFSs would facilitate dispersal over selection, thereby preventing the establishment of distinct benthic biofilm communities.


We sampled benthic biofilms and streamwater from 54 GFSs across the European Alps, the Scandinavian Mountains, and the Himalayas. Studying GFSs from the three mountain ranges allowed us to draw more generalizable conclusions on their microbial communities and assembly processes. Within each GFS, we sampled benthic sandy sediments from three independent patches downstream from the glacier snout. In addition, water samples were collected from each GFS and sequentially filtered in the field to isolate different fractions of the suspended bacterial assemblage (see Materials and Methods). Thereby, while we were able to sample UC alone, AC could not be separated from UC; we therefore treated AC along with UC as a composite sample (UC+AC). We considered UC and UC+AC together as the assemblages in the streamwater from which benthic biofilms may assemble.

Taxonomic composition.

We analyzed the taxonomic composition of benthic biofilms and water assemblages across all three mountain ranges. We found that a large number of ASVs were specific to the benthic biofilms (31,251 ASVs), UC (25,086 ASVs), and UC+AC (12,358 ASVs) (Fig. 1A). However, when taking into account the relative abundance of ASVs, we noted that only a small proportion of taxa was specific to the benthic biofilms (4.7%) and the source assemblages (5.5% and 2% for the UC and UC+AC, respectively). This suggests that a low number of highly abundant taxa was shared across the benthic biofilms and the assemblages suspended in the streamwater (i.e., 6,878 taxa accounting for 68% of the ASVs), while an elevated number of rare taxa (i.e., low ratio of proportion to relative abundance) was specific to each.
FIG 1 Taxonomic composition of benthic biofilms and source assemblages. (A) The upset plot depicts the number of unique and shared ASVs among benthic biofilms, unattached (UC), and attached cells (UC+AC) from the overlying streamwater and the relative abundance characterizing each set. (B) The stacked bar plot illustrates the average relative abundance of shared and/or unique ASVs across benthic biofilms and water assemblages. To increase clarity on the barplot, features were agglomerated at the family level.
For instance, only three phyla dominated the bacterial assemblages across the entire data set, including Proteobacteria (49.6% ± 14.2%), Patescibacteria (11.6 ± 15.8), and Bacteroidota (9.9% ± 6.1%; Table S1). Benthic biofilms were dominated by members from the families Comamonadaceae (17% ± 7.2%), Sphingomonadaceae (8.3% ± 4.9%), Methylophilaceae (4.3% ± 4.9%), and Chitinophagaceae (3.3% ± 2%), while members of Candidatus Nomurabacteria (11.2% ± 5.7%), Comamonadaceae (9.7% ± 8.9%) and Sporichthyaceae (4.8% ± 4.8%) dominated UC, and lineages from the families Comamonadaceae (15% ± 7.6%), Candidatus Nomurabacteria (4.5% ± 2.8%), Flavobacteriaceae (4.1% ± 4.4%), and Methylophilaceae (3.4% ± 4.3%) mainly populated UC+AC water assemblages (Fig. 1B; Tables S2 to 4).

Diversity and compositional patterns.

First, we tested for differences in alpha diversity between benthic biofilms (from three patches in each GFS) and their potential streamwater source assemblages (UC+AC, UC). We found that amplicon sequence variant (ASV) richness was significantly lower (analysis of variance [ANOVA]; Ffraction = 6.4; pfraction = 0.002; Tukey HSD, ppairwise = 0.002) in the benthic biofilms (1497 ± 425) compared with the UC fraction (1775 ± 634) across the three mountain ranges (Fig. 2A; Table S5), but not than those of UC+AC (Tukey HSD, ppairwise = 0.22). Although not statistically significant, the average ASV richness was higher in UC than in UC+AC, which could be attributed to different filtration yields of sample volumes and hence cell numbers (see Materials and Methods). Similar differences in richness between streamwater assemblages and biofilms growing on pebbles (8) and rocks (23) have been previously reported from other GFSs and low-land rivers (9). Collectively, these patterns support the idea that the elevated richness in the streamwater is attributable to the fact that GFSs “collect” bacteria from various catchment sources (e.g., glacier, soils, rock, snow).
FIG 2 Patterns of alpha and beta diversity for benthic biofilms and water assemblages. (A) Richness of ASV for benthic biofilms (n = 142), composite samples of unattached cells (UC) and cells attached to suspended particles (AC) suspended in the streamwater (AC+UC; n = 46), and UC (n = 45) alone. Box plots show the median (horizontal line), interquartile range (box height), 1.5x beyond the interquartile range (whiskers), and outliers. Letters indicate significant differences based on pairwise comparisons (P < 0.05). (B) Differences in community composition shown in the nonmetric multidimensional scaling (NMDS) space as based on Bray-Curtis dissimilarity.
Next, we explored patterns of beta diversity to assess possible compositional differences between the benthic biofilms and the streamwater assemblages. Using nonmetric multidimensional scaling (NMDS) analysis based on Bray-Curtis dissimilarity, we found that their composition differed significantly (anova.manyglm, Dev = 41923, pfraction x mountain range = 0.01; ppairwise - benthic biofilm x water assemblages < 0.03 for all tests) and consistently across all GFSs from the three mountain ranges (Fig. 2B; Table S6, S7). Interestingly, the composition of benthic biofilm communities was more similar to the composition of UC+AC (average BCbenthic biofilm versus UC+AC = 0.91) than UC alone (average BCbenthic biofilm versus UC = 0.96). We attribute this pattern to the AC assemblage potentially containing attached bacteria that are more likely to form suspended biofilms compared to UC. Furthermore, AC may be transported more easily from the streamwater toward the benthic biofilms because of the settling properties of the sediment particles they are attached to. Differing patterns in community composition between streamwater and biofilms collected from pebbles (8) and rocks (23) were observed in GFSs and nonglacier headwater streams (9). In line with its elevated alpha diversity, UC also exhibited the highest level of compositional variability compared to benthic biofilms (betadisper, F = 4.9, pfraction = 0.01; Tukey HSD, P = 0.007; Table S8, 9), which further corroborates the notion of various sources contributing to this fraction and differing among GFSs. That said, compositional differences were mainly attributable to changes in the relative abundance of the most abundant ASVs and presence/absences of rare ASVs that characterize benthic biofilms, UC+AC, and UC (Fig. 1; Table S1 to 4). A similar role of abundant and rare taxa for differences in community composition between benthic and suspended bacterial assemblages was reported from a low-land stream (9).

Differential abundance.

The differences in diversity and community composition between the benthic biofilms and suspended source assemblages prompted us to further investigate their taxonomic composition. Using differential abundance analyses (Fig. 3; Table S10), we found bacterial families differentially enriched in benthic biofilms compared to UC+AC (Log2 median ratio benthic biofilm versus UC+AC < –0.66 and padj < 0.05) and UC (Log2 median ratio benthic biofilm versus UC < –1.3, padj < 0.04), respectively. It is worth noting that most families occurring in greater abundance in the benthic biofilms were also identified in UC+AC and UC, albeit at lower abundances.
FIG 3 Differential abundance analyses highlighting taxonomic differences between benthic biofilms and streamwater assemblages across the bacterial tree. (A) Taxonomic tree showing bacterial taxa that are on average over-represented in the benthic biofilms compared with the unattached cells (UC) in the streamwater. (B) Taxonomic tree showing bacterial taxa that are on average over-represented in the benthic biofilms compared with the unattached and attached cells (UC+AC) in the streamwater. Size and color of the edges and nodes are directly correlated with the abundance of bacteria within each fraction.
As shown on Fig. 3, we found 45 families to be enriched in the benthic biofilms compared with both UC+AC and UC (Table S10) spanning 12 phyla commonly reported from soils and aquatic ecosystems, including various glacial environments (2729). For instance, Caulobacteraceae, Hyphomicrobiaceae, Rhodobacteraceae, Sphingomonadaceae, and Nitrosomonadaceae belonging to Proteobacteria prevailed at moderate to high relative abundance in the benthic biofilm (0.09% up to 8.3%; Table S2). Some of these taxa are known to produce exo-lipopolysaccharides (30), which are important components of the biofilm extracellular matrix (31). Along with Proteobacteria, we found families that were moderately overrepresented (0.7% to 1.4%) in the benthic biofilms compared with UC+AC and UC and affiliated with the phyla Verrucomicrobia, Nitrospirota, Bacteroidota, and Planctomycetes (Table S10). We also detected 19 families with low abundances (<1%; Table S2) that were overrepresented in the benthic biofilms and belonging to the Acidobacteriota and Chloroflexi phyla (Table S10). The phylum Acidobacteriota contains psychrotolerant, mesophilic and slow-growing taxa, primarily in nutrient-poor soils (32), while Chloroflexi encompasses filamentous photoheterotrophic bacteria known to degrade complex organic compounds (33, 34).
Comparing benthic biofilms with UC only, we found 58 families specifically overrepresented in the former (Fig. 3A; Table S10), some of them reported from soils and various glacial environments, and involved in biofilm formation and putatively important biogeochemical processes in streams (e.g., sulfur, nitrogen, carbon cycling) (8, 9, 19, 3537). These families included members of Proteobacteria (e.g., Comamonadaceae, Xanthomonadaceae), Bacteroidota (e.g., Flavobacteriaceae, Chitinophagaceae), and Myxococcota (e.g., Myxococcaceae). A few of these families, such as Comamonadaceae and Chitinophagaceae, were among the most abundant in the benthic biofilms (e.g., 17% for Comamonadaceae).
Taxa that were over-represented in the UC+AC and UC assemblages compared with the benthic biofilms included numerous members from the Patescibacteria super phylum, but also Proteobacteria (e.g., Holosporaceae, Rickettsiaceae, Pseudomonadaceae), Chlamydiota, Actinobacteria, and families from the order Bacillales, for instance. Some of these taxa are potential spore formers (38), others harbor ultra-small cells and may associate (e.g., episymbiosis) with other microorganisms for nutrient supply (39) (Table S10).

Community assembly processes.

Our analyses on diversity, composition, and taxonomic patterns indicate that benthic biofilms are distinct compared with the overlying source assemblages in the streamwater. To further substantiate this notion and assess potential processes (i.e., stochastic versus deterministic) driving the assembly of the benthic biofilms, we applied a neutral model (24) and a phylogenetic/compositional turnover framework (25, 26) (see Materials and Methods). Our null hypothesis was that there are no ecological differences between the benthic biofilm and source assemblages, such that the observed differences among them could be explained by source-sink dynamics and dispersal.
The Sloan’s neutral community model that we first applied (26) assumes that there are no constraints in dispersal across a given metacommunity, such that the relationship between the relative abundance of a given taxon and its occurrence frequency is described by a beta distribution. This model can be also applied by considering a taxon’s abundance in a source metacommunity and its occurrence frequency in the target metacommunity, respectively, to test for source-sink hypotheses. Results from the neutral model showed that, across the three mountain ranges, the composition of benthic biofilms is not driven by a source-sink relationship with UC+AC or UC (Fig.S1, 2; Table S11).
The phylogenetic/compositional turnover framework that we next applied is based on null modeling, and assigns the patterns of phylogenetic turnover that are distinct from the null distribution to niche-related processes and the respective patterns of compositional turnover to dispersal-related or stochastic processes (25, 26). Results from this framework revealed lower than expected phylogenetic turnover between benthic biofilms and UC+AC across all three mountain ranges (Fig. 4A). This highlights homogenous selection as the dominant assembly process of the biofilm communities because lower than expected phylogenetic turnover suggests higher than expected ecological similarities among biofilm communities (under the notion that phylogenetic similarity correlates to ecological similarity for bacteria) (40). This analysis also revealed a number of processes, mainly variable selection (higher than expected phylogenetic turnover) and dispersal limitation (higher than expected compositional turnover), to drive differences in community assembly between benthic biofilms and UC (Fig. 4B). The lack of dominance of homogenizing dispersal (lower than expected compositional turnover) is also evident in both the comparisons, substantiating the results obtained from Sloan’s neutral community model.
FIG 4 Dominant assembly processes regarding sample pairs of benthic biofilms compared with the assemblages in the streamwater for glacier-fed streams from each mountain range. (A) Results regarding comparisons of benthic biofilms versus unattached and attached cells (UC+AC). (B) Results regarding comparisons of benthic biofilms versus unattached cells (UC) only.
Collectively, our results above reject our initial null hypothesis; dispersal from UC (partially also from UC+AC) is not effectively homogenizing community structure of the benthic biofilms or imposing source-sink dynamics. These findings hint toward strong species sorting and competitive exclusion in benthic biofilms as also inferred by previous studies from various stream ecosystems (8, 9, 41). This might be associated with the fact that, within benthic biofilms, the residence time for bacterial cells is higher than in the overlying streamwater, which enables them to produce extracellular polymeric substances, construct niches, and engage in interactions with other taxa, even across kingdoms (31, 42). The biofilm mode of life therefore involves niche-related processes, which would also explain the differences in alpha and beta diversity in GFS benthic biofilms and their source assemblages in the streamwater. In contrast, the more-stochastic assembly in the UC fraction is consistent with findings from other works showing that stochasticity can sometimes lead to an increase in alpha diversity (5, 43), and agrees with the notion that the GFS streamwater probably “collects” microorganisms from various catchment sources.


Our results show that benthic biofilms across a wide range of GFSs differ in alpha and beta diversity from the assemblages in the overlying streamwater. This is striking given the continuous vertical mixing of suspended cells (attached to particles and unattached) with the benthic sediments in these low-submergence and highly turbulent streams. Our results revealed deterministic processes favoring the presence of specific taxa and driving community assembly of the benthic biofilms. Collectively, these findings highlight the benthic sediments of GFSs as a distinct habitat with specific biofilm communities and reveal a “vertical” dimension of metacommunity dynamics in GFSs. These findings are critical given that GFS ecosystems are rapidly changing owing to climate-induced glacier shrinkage (22, 44). As glaciers shrink, the sedimentary environment in GFSs is expected to change as well, and so will the sources of microorganisms to the GFSs. Further large-scale studies are required to better understand and predict these climate-induced impacts on the microbial life in GFSs.


Study sites and sample collection.

Benthic biofilm and water samples were collected from 54 GFSs across the European Alps, the Scandinavian Mountains, and the Himalayas (Table S12; NAlps = 27, NScandinavian Mountains = 10, N Himalayas = 17) between August 2020 and May 2021. The glacier catchment area varied from 0.4 to 105 km2. Within each GFS, we targeted three independent patches located within a 5-m radius, as close as possible to the glacier snout (1 m to 910 m). From each patch, we sampled up to 30 g of benthic sediments (per sample, with one sample per patch) by collecting and sieving (3.15 mm to 250 μm, Retsch, Woven Wire Mesh Sieve) the upper 5 cm of the stream bed with flame-sterilized metal tools. In addition, water samples were collected from each GFS to compare suspended bacterial assemblages with the benthic biofilms. Water samples were filtered to isolate two different size fractions (one sample per fraction for each GFS), aimed at isolating different microbial assemblages within the water column. Specifically, the UC fraction was isolated through a stacked filtration system (1 μm Polycap cartridge fitted onto a 0.2 μm Sterivex) by filtering on average 1,210 ± 745 mL. The corresponding UC+AC was only filtered through a 0.2 μm Sterivex filter with an average volume of 413 ± 361 mL and isolated both AC and UC. From microscopic analyses of raw unfiltered water from glacier samples, we noted that bacterial cells showed an average length of 0.65 ± 0.5 μm (Fig. S3). Samples were immediately frozen in liquid nitrogen and subsequently stored in –80°C prior to and following shipping.

DNA extraction and metabarcoding library preparation.

DNA extraction was performed on sediments using a phenol-chloroform protocol adjusted to the specificities of our samples (45). Regarding DNA extraction from water filters, we used a phenol-chloroform modified version of the protocol by Massana et al. (46). DNA quantification was performed for all samples with the Qubit dsDNA HS kit (Invitrogen). As numerous sediment samples exhibited low DNA yields, all samples were diluted to a final concentration of ≤ 2 to 3 ng/μL to avoid PCR bias. The V3-V4 hypervariable regions of the bacterial 16S rRNA gene were amplified using primers 341f (5′-CCTACGGGNGGCWGCAG-3′) and 785r (5′-GACTACHVGGGTATCTAATCC-3′), synthesized with the appropriate overhang adapter sequences for the Illumina MiSeq system. The KAPA HiFi DNA polymerase (Hot Start and Ready Mix formulation) was used in a 25-μL amplification reaction containing 1 × PCR buffer, 1 μM each primer, 0.48 μg/μL BSA, and 1.0 μL of template DNA. Amplification was performed in a Biometra Trio (Biometra) instrument. The thermal conditions applied after an initial denaturation at 95°C for 3 min were 94°C for 30 s, 55°C for 30 s, and 72°C for 30 s for 25 cycles followed by a final extension at 72°C for 5 min. Amplification was verified on a 1.5% agarose gel. Amplicon libraries were prepared according to the MiSeq manufacturer’s protocol. In brief, a second PCR was performed for the addition of dual indices to the purified amplicon PCR products. This allowed extensive multiplexing of samples on a single sequencing lane of the MiSeq (Illumina) platform after quantification and normalization. Samples were sequenced using a 300 base paired-end protocol in the Biological Core Lab of Red Sea Research Center of King Abdullah University of Science and Technology in Saudi Arabia.

Amplicon sequence data processing.

A total of 234 amplicon sequence libraries were generated, including sediment and filters, after 26 samples were lost following DNA extraction due to low biomass and an additional 10 samples lost during transportation in the Himalayas. Paired-end sequencing generated a total of 42,896,962 reads across 234 samples, including an average read number of 170,450 and 158,717 for sediment and water samples, respectively. Raw sequences were trimmed using the plugin Trimmomatic v.0.36 as outlined in Fodelianakis et al. (41). Briefly, reads were truncated in 4 bp sliding windows when mean quality dropped below a Phred score of 15. We removed the three leading and trailing nucleotides and discarded reads shorter than 200 bp. All subsequent sequence analyses were processed using the Quantitative Insights into Microbial Ecology 2 (QIIME2 v.2020.8) (47) workflow. We used the demux plugin to visualize interactive quality plots and assess read quality. Using DADA2 (48) with default parameters, primers were trimmed, and reads were denoised and joined to produce ASVs. Specifically, 17 and 21 nucleotides were removed at the beginning of each forward and reverse reads. Forward and reverse reads were truncated at 280 and 220 bases, respectively. These different steps resulted in a total of 27,682,066 reads across 234 samples. We then assigned taxonomy against the SILVA reference database (v.138) (49) using the algorithm classify-sklearn (50). A phylogenetic tree was computed with the alignment and phylogeny plugins for further downstream analyses. The biom table including the taxonomic counts, the phylogenetic tree and related metadata were imported in R for further statistical analyses. Mitochondria, chloroplasts, Archaea related reads, sequences related to bacteria that did not map any classified phyla, together with singletons and doubletons were removed, further resulting in a total of 25,781,108 reads across 234 samples. Rarefaction curves were generated (Fig. S4). We used a random resampling method without replacement to rarefy our ASV table. Rarefaction threshold was set at 17,309 reads per sample across 233 samples, after one water sample was discarded due to low sampling depth (<1,000 reads).

Data analyses.

All statistical analyses were performed in R (4.1.0). First, we investigated patterns of alpha diversity across the data set using the ASV richness computed with the function plot_richness within the phyloseq package (v.1.27.6) (51). We assessed the effects of the “fraction,” “mountain range,” and their interaction on ASV richness using a two-way ANOVA and pairwise comparisons among group levels were further tested using Tukey HSD present in the package stats (4.1.0). Data residuals were tested for normality and homoscedasticity using Shapiro-Wilk and Levene tests.
To display the number and proportion of ASVs shared between fractions, we computed an UpSet plot using the function upset in the package UpSetR (v.1.4.0) (52). To ensure a clear representation of the stacked barplot, taxa were agglomerated at the family level using the function tax_glom (including the command NArm = F) within the phyloseq package. To illustrate potential shifts in bacterial assemblages as a function of “fraction” and “mountain range,” we computed a nonmetric multidimensional scaling (NMDS) on the Bray-Curtis dissimilarity matrix (including 999 permutations) using the function plot_ordination in phyloseq. We assessed differences in bacterial community composition as a function of the “fraction,” “mountain range,” and their interaction using a negative binomial multivariate generalized linear model (functions manyglm and anova.manyglm in the package mvabund [v.4.2.1]) (53). To account for multiple comparisons, pairwise differences were tested and P-values were adjusted using a free stepdown resampling procedure.
We tested the effects of “fraction” and “mountain range” on bacterial compositional variability using an analysis of multivariate homogeneity of group dispersions as a function of “fraction” and “mountain range” on the Bray-Curtis dissimilarity matrices using the function betadisper in vegan (54). When significant, Tukey HSD tests were used to test for significance of pairwise comparisons.
The package Metacoder (v.0.3.5) (55) was used to compute differential abundance analyses between the different fractions and to plot the related taxonomic data through a differential heat tree. The unrarefied data set was loaded prior to running the Metacoder pipeline and low-abundance counts were discarded (i.e., min_count was set to 5). Wilcoxon Rank Sum tests were computed to test for differences between the median abundance of samples across fractions. P-values were adjusted using the Holm method (56).
To investigate source-sink relationships between benthic biofilms and water fractions, we utilized the Sloan’s neutral community model (24) using the R code developed in Burns et al. (57). Finally, to quantify the dominant processes driving community assembly, we first computed the phylogenetic turnover (β-nearest taxon index: βNTI) and then assessed the compositional turnover (Raup-Crick distances based on the Bray-Curtis dissimilarity [RCBC]) at the community level following the framework by Stegen et al. (25, 26). Briefly, through this framework, differences between two communities (i.e., benthic biofilms and UC or UC+AC) are either assigned to selection (homogenous or heterogenous), dispersal (homogenizing or limiting), or to a lack of any dominant process. To evaluate the influence of selection, we first examined phylogenetic turnover at the community level. Under the notion that phylogenetic similarity correlates with ecological similarity, as is commonly the case for bacterial communities (40), lower than expected phylogenetic turnover between two given communities suggests higher than expected ecological similarities and thus similar niche-related processes acting in both communities (termed “homogenous selection”) and vice versa. We determined the community-wise phylogenetic turnover between two given communities via the z-score (e.g., β-nearest taxon index or βNTI) of the β-mean nearest taxon distance from a null distribution of the same metric. Values of βNTI less than –2 suggest that homogenous selection between the compared communities causes higher-than-expected phylogenetic similarity. In contrast, βNTI scores greater than +2 are an indication of heterogeneous selection. Community pairs with |βNTI| less than 2 (nonsignificantly different phylogenetic turnover to that expected by chance) are then further investigated in terms of compositional turnover. Here, the notion is that niche-unrelated processes “blind to phylogeny” such as (passive) dispersal via some abiotic vector (in our case water), could still homogenize two given communities if they are strong enough, or cause compositional differences if they are weak enough (limiting). To examine such patterns, we utilized the RCBC metric. Values of RCBC less than –0.95 are attributed to strong dispersal (termed “homogenizing dispersal”) while RCBC values greater than 0.95 are indicative of dispersal limitation.
All data presented in the Results and Discussion section are described as mean ± standard deviation.

Data availability.

Raw sequence data have been deposited in the Sequence Read Archives (SUB10671583, SUB11121237, SUB10689269) with NCBI BioProject accession no. PRJNA781406. The metadata relative to GFS sediment and water samples are available on the NOMIS database platform: nomis-data.epfl.ch and R code underlying statistical analyses and models has been deposited to github.com/laylaeb/RIVER_GFS_WS. Additional information related to the NOMIS Foundation project “Vanishing Glaciers” could be found on the official website: https://www.glacierstreams.ch/.


We would like to acknowledge the support from Subodh Sharma and Prayon Joshi from Kathmandu University (Nepal), Antoine Rabat, Jean Baptiste Bosson, Brad Carlson and Richard Bonet for providing help when conducting sampling in the European Alps as well as Andreas Schomacker, Steffen Leth Jørgensen, Ida Marie Gabrielsen, Jacob Clement Yde when sampling in the Scandinavian Mountains.
This work was supported by The NOMIS Foundation project “Vanishing Glaciers” to T.J.B. S.B.B. was supported by the Synergia grant (CRSII5_180241: Swiss National Science Foundation) to T.J.B. D.D. acknowledges the financial support of King Abdullah University and Technology (KAUST) through the baseline research fund. We also thank two anonymous reviewers whose comments greatly improved the manuscript.
We declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
L.E. and T.J.B. conceived the project and designed the experiment; Mi.S., V.D.S., Ma.S., M.T., T.J.K., and T.J.B. conducted field work; P.P. and E.O. conducted DNA extraction and library preparation, and R.M. and D.D. performed the sequencing; H.P. and V.T. conducted microscopic analyses; L.E. and S.F. analyzed the data; L.E., T.J.K., S.F., and T.J.B. wrote the manuscript with inputs and editing from all coauthors. All authors contributed to the article and approved the submitted version.

Supplemental Material

File (aem.00421-22-s0001.xlsx)
File (aem.00421-22-s0002.pdf)
ASM does not own the copyrights to Supplemental Material that may be linked to, or accessed through, an article. The authors have granted ASM a non-exclusive, world-wide license to publish the Supplemental Material files. Please contact the corresponding author directly for reuse.


Rinaldo A, Gatto M, Rodríguez-Iturbe I. 2020. River networks as ecological corridors: species, populations, pathogens. Cambridge University Press.
Ruiz-González C, Niño-García JP, del Giorgio PA. 2015. Terrestrial origin of bacterial communities in complex boreal freshwater networks. Ecol Lett 18:1198–1206.
Savio D, Sinclair L, Ijaz UZ, Parajka J, Reischer GH, Stadler P, Blaschke AP, Blöschl G, Mach RL, Kirschner AKT, Farnleitner AH, Eiler A. 2015. Bacterial diversity along a 2600 km river continuum. Environ Microbiol 17:4994–5007.
Crump BC, Amaral-Zettler LA, Kling GW. 2012. Microbial diversity in arctic freshwaters is structured by inoculation of microbes from soils. ISME J 6:1629–1639.
Liu K, Liu Y, Hu A, Wang F, Zhang Z, Yan Q, Ji M, Vick‐Majors TJ. 2021. Fate of glacier surface snow‐originating bacteria in the glacier‐fed hydrologic continuums. Environ Microbiol 23:6450–6462.
Heino J, Melo AS, Siqueira T, Soininen J, Valanko S, Bini LM. 2015. Metacommunity organisation, spatial extent and dispersal in aquatic systems: patterns, processes and prospects. Freshw Biol 60:845–869.
Gu Z, Liu K, Pedersen MW, Wang F, Chen Y, Zeng C, Liu Y. 2021. Community assembly processes underlying the temporal dynamics of glacial stream and lake bacterial communities. Sci Total Environ 761:143178.
Wilhelm L, Singer GA, Fasching C, Battin TJ, Besemer K. 2013. Microbial biodiversity in glacier-fed streams. ISME J 7:1651–1660.
Besemer K, Peter H, Logue JB, Langenheder S, Lindström ES, Tranvik LJ, Battin TJ. 2012. Unraveling assembly of stream biofilm communities. ISME J 6:1459–1468.
Anderson SP. 2007. Biogeochemistry of glacial landscape systems. Annu Rev Earth Planet Sci 35:375–399.
Milner AM, Brittain JE, Brown LE, Hannah DM. 2010. Water sources and habitat of alpine streams, p. 175–191. Alpine waters. Springer, Heidelberg, Berlin.
Milner AM, Petts GE. 1994. Glacial rivers: physical habitat and ecology. Freshwater Biol 32:295–307.
Jacobsen D, Dangles O. 2012. Environmental harshness and global richness patterns in glacier‐fed streams. Glob Ecol Biogeogr 21:647–656.
Hotaling S, Hood E, Hamilton TL. 2017. Microbial ecology of mountain glacier ecosystems: biodiversity, ecological connections and implications of a warming climate. Environ Microbiol 19:2935–2948.
Cauvy-Fraunié S, Dangles O. 2019. A global synthesis of biodiversity responses to glacier retreat. Nat Ecol Evol 3:1675–1685.
Van Horn DJ, Wolf CR, Colman DR, Jiang X, Kohler TJ, McKnight DM, Stanish LF, Yazzie T, Takacs-Vesbach CD. 2016. Patterns of bacterial biodiversity in the glacial meltwater streams of the McMurdo Dry Valleys, Antarctica. FEMS Microbiol Ecol 92:fiw148.
Battin TJ, Kaplan LA, Newbold JD, Cheng X, Hansen C. 2003. Effects of current velocity on the nascent architecture of stream microbial biofilms. Appl Environ Microbiol 69:5443–5452.
Busi SB, Bourquin M, Fodelianakis S, Michoud G, Kohler TJ, Peter H, Pramateftaki P, Styllas M, Tolosano M, De Staercke V. 2021. Genomic and metabolic adaptations of biofilms to ecological windows of opportunities in glacier-fed streams.
Wilhelm L, Besemer K, Fasching C, Urich T, Singer GA, Quince C, Battin TJ. 2014. Rare but active taxa contribute to community dynamics of benthic biofilms in glacier‐fed streams. Environ Microbiol 16:2514–2524.
Boix Canadell M, Gómez‐Gener L, Ulseth AJ, Clémençon M, Lane SN, Battin TJ. 2021. Regimes of primary production and their drivers in Alpine streams. Freshw Biol 66:1449–1463.
Leibold MA, Holyoak M, Mouquet N, Amarasekare P, Chase JM, Hoopes MF, Holt RD, Shurin JB, Law R, Tilman D, Loreau M, Gonzalez A. 2004. The metacommunity concept: a framework for multi‐scale community ecology. Ecol Lett 7:601–613.
Zemp M, Frey H, Gärtner-Roer I, Nussbaumer SU, Hoelzle M, Paul F, Haeberli W, Denzinger F, Ahlstrøm AP, Anderson B, Bajracharya S, Baroni C, Braun LN, Cáceres BE, Casassa G, Cobos G, Dávila LR, Delgado Granados H, Demuth MN, Espizua L, Fischer A, Fujita K, Gadek B, Ghazanfar A, Ove Hagen J, Holmlund P, Karimi N, Li Z, Pelto M, Pitte P, Popovnin VV, Portocarrero CA, Prinz R, Sangewar CV, Severskiy I, Sigurđsson O, Soruco A, Usubaliev R, Vincent C. 2015. Historically unprecedented global glacier decline in the early 21st century. J Glaciol 61:745–762.
Hotaling S, Foley ME, Zeglin LH, Finn DS, Tronstad LM, Giersch JJ, Muhlfeld CC, Weisrock DW. 2019. Microbial assemblages reflect environmental heterogeneity in alpine streams. Glob Chang Biol 25:2576–2590.
Sloan WT, Lunn M, Woodcock S, Head IM, Nee S, Curtis TP. 2006. Quantifying the roles of immigration and chance in shaping prokaryote community structure. Environ Microbiol 8:732–740.
Stegen JC, Lin X, Fredrickson JK, Konopka AE. 2015. Estimating and mapping ecological processes influencing microbial community assembly. Front Microbiol 6:370.
Stegen JC, Lin X, Fredrickson JK, Chen X, Kennedy DW, Murray CJ, Rockhold ML, Konopka A. 2013. Quantifying community assembly processes and identifying features that impose them. ISME J 7:2069–2079.
Tolotti M, Cerasino L, Donati C, Pindo M, Rogora M, Seppi R, Albanese D. 2020. Alpine headwaters emerging from glaciers and rock glaciers host different bacterial communities: ecological implications for the future. Sci Total Environ 717:137101.
Krauze P, Wagner D, Yang S, Spinola D, Kühn P. 2021. Influence of prokaryotic microorganisms on initial soil formation along a glacier forefield on King George Island, maritime Antarctica. Sci Rep 11:13135.
Ren Z, Gao H, Elser JJ. 2017. Longitudinal variation of microbial communities in benthic biofilms and association with hydrological and physicochemical conditions in glacier-fed streams. Freshw Sci 36:479–490.
Vuko M, Cania B, Vogel C, Kublik S, Schloter M, Schulz S. 2020. Shifts in reclamation management strategies shape the role of exopolysaccharide and lipopolysaccharide‐producing bacteria during soil formation. Microb Biotechnol 13:584–598.
Flemming H-C, Wingender J. 2010. The biofilm matrix. Nat Rev Microbiol 8:623–633.
Dedysh SN, Damsté JSS. 2018. Acidobacteria. eLS, p. 1–10. John Wiley & Sons, Ltd., Chichester, UK.
Liu Y, Vick-Majors TJ, Priscu JC, Yao T, Kang S, Liu K, Cong Z, Xiong J, Li Y. 2017. Biogeography of cryoconite bacterial communities on glaciers of the Tibetan Plateau. FEMS Microbiol Ecol 93:fix072.
Speirs LBM, Rice DTF, Petrovski S, Seviour RJ. 2019. The phylogeny, biodiversity, and ecology of the chloroflexi in activated sludge. Front Microbiol 10.
Battin TJ, Besemer K, Bengtsson MM, Romani AM, Packmann AI. 2016. The ecology and biogeochemistry of stream biofilms. Nat Rev Microbiol 14:251–263.
Nemergut DR, Anderson SP, Cleveland CC, Martin AP, Miller AE, Seimon A, Schmidt SK. 2007. Microbial community succession in an unvegetated, recently deglaciated soil. Microb Ecol 53:110–122.
Zhang L, Delgado-Baquerizo M, Shi Y, Liu X, Yang Y, Chu H. 2021. Co-existing water and sediment bacteria are driven by contrasting environmental factors across glacier-fed aquatic systems. Water Res 198:117139.
Nicholson WL, Munakata N, Horneck G, Melosh HJ, Setlow P. 2000. Resistance of Bacillus endospores to extreme terrestrial and extraterrestrial environments. Microbiol Mol Biol Rev 64:548–572.
Tian R, Ning D, He Z, Zhang P, Spencer SJ, Gao S, Shi W, Wu L, Zhang Y, Yang Y, Adams BG, Rocha AM, Detienne BL, Lowe KA, Joyner DC, Klingeman DM, Arkin AP, Fields MW, Hazen TC, Stahl DA, Alm EJ, Zhou J. 2020. Small and mighty: adaptation of superphylum Patescibacteria to groundwater environment drives their genome simplicity. Microbiome 8:1–15.
Dini-Andreote F, Stegen JC, Van Elsas JD, Salles JF. 2015. Disentangling mechanisms that mediate the balance between stochastic and deterministic processes in microbial succession. Proc Natl Acad Sci USA 112:E1326–E1332.
Fodelianakis S, Washburne AD, Bourquin M, Pramateftaki P, Kohler TJ, Styllas M, Tolosano M, De Staercke V, Schön M, Busi SB, Brandani J, Wilmes P, Peter H, Battin TJ. 2022. Microdiversity characterizes prevalent phylogenetic clades in the glacier-fed stream microbiome. ISME J 16:666–610.
Hall-Stoodley L, Costerton JW, Stoodley P. 2004. Bacterial biofilms: from the natural environment to infectious diseases. Nat Rev Microbiol 2:95–108.
Xun W, Li W, Xiong W, Ren Y, Liu Y, Miao Y, Xu Z, Zhang N, Shen Q, Zhang R. 2019. Diversity-triggered deterministic bacterial assembly constrains community functions. Nat Commun 10:1–10.
Zemp M, Haeberli W. 2007. Glaciers and ice caps. Part I: global overview and outlook. Part II: glacier changes around the world. UNEP Glob Outlook Ice Snow :115–152.
Busi SB, Pramateftaki P, Brandani J, Fodelianakis S, Peter H, Halder R, Wilmes P, Battin TJ. 2020. Optimised biomolecular extraction for metagenomic analysis of microbial biofilms from high-mountain streams. PeerJ 8:e9973.
Massana R, Murray AE, Preston CM, DeLong EF. 1997. Vertical distribution and phylogenetic characterization of marine planktonic Archaea in the Santa Barbara Channel. Appl Environ Microbiol 63:50–56.
Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet C, Al-Ghalith GA, Alexander H, Alm EJ, Arumugam M, Asnicar F. 2018. QIIME 2: reproducible, interactive, scalable, and extensible microbiome data science. PeerJ Preprints.
Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. 2016. DADA2: high-resolution sample inference from Illumina amplicon data. Nat Methods 13:581–583.
Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, Peplies J, Glöckner FO. 2013. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res 41:D590–D596.
Bokulich NA, Kaehler BD, Rideout JR, Dillon M, Bolyen E, Knight R, Huttley GA, Caporaso JG. 2018. Optimizing taxonomic classification of marker-gene amplicon sequences with QIIME 2’s q2-feature-classifier plugin. Microbiome 6:1–17.
McMurdie PJ, Holmes S. 2013. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS One 8:e61217.
Conway JR, Lex A, Gehlenborg N. 2017. UpSetR: an R package for the visualization of intersecting sets and their properties. Bioinformatics 33:2938–2940.
Wang YI, Naumann U, Wright ST, Warton DI. 2012. mvabund–an R package for model-based analysis of multivariate abundance data. Methods Ecol Evol 3:471–474.
Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, McGlinn D, Minchin PR, O’Hara RB, Simpson GL, Solymos P. 2020. vegan: community ecology package. R package version 2.5–6.
Foster ZSL, Sharpton TJ, Grünwald NJ. 2017. Metacoder: an R package for visualization and manipulation of community taxonomic diversity data. PLoS Comput Biol 13:e1005404.
Holm S. 1979. A simple sequentially rejective multiple test procedure. Scand J Stat 6:65–70.
Burns AR, Stephens WZ, Stagaman K, Wong S, Rawls JF, Guillemin K, Bohannan BJM. 2016. Contribution of neutral processes to the assembly of gut microbial communities in the zebrafish over host development. ISME J 10:655–664.

Information & Contributors


Published In

cover image Applied and Environmental Microbiology
Applied and Environmental Microbiology
Volume 88Number 1228 June 2022
eLocator: e00421-22
Editor: Knut Rudi, Norwegian University of Life Sciences
PubMed: 35674429


Received: 8 March 2022
Accepted: 17 May 2022
Published online: 8 June 2022


Request permissions for this article.


  1. benthic biofilms
  2. community composition
  3. community assembly
  4. glacier-fed streams



River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland
Stilianos Fodelianakis
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland
Tyler J. Kohler
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland
Massimo Bourquin
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland
Jade Brandani
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland
Susheel Bhanu Busi
Systems Ecology group, Luxembourg Centre for Systems Biomedicine, University of Luxembourg, Esch-sur-Alzette, Luxembourg
Daniele Daffonchio
Biological and Environmental Sciences and Engineering Division (BESE), King Abdullah University of Science and Technology (KAUST), Thuwal, Saudi Arabia
Vincent De Staercke
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland
Biological and Environmental Sciences and Engineering Division (BESE), King Abdullah University of Science and Technology (KAUST), Thuwal, Saudi Arabia
Grégoire Michoud
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland
Emmy Oppliger
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland
Paraskevi Pramateftaki
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland
Martina Schön
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland
Michail Styllas
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland
Virginia Tadei
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland
Matteo Tolosano
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland
Tom J. Battin [email protected]
River Ecosystems Laboratory, Alpine and Polar Environmental Research Center, Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland


Knut Rudi
Norwegian University of Life Sciences


The authors declare no conflict of interest.

Metrics & Citations



  • For recently published articles, the TOTAL download count will appear as zero until a new month starts.
  • There is a 3- to 4-day delay in article usage, so article usage will not appear immediately after publication.
  • Citation counts come from the Crossref Cited by service.


If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. For an editable text file, please select Medlars format which will download as a .txt file. Simply select your manager software from the list below and click Download.

View Options

Figures and Media






Share the article link

Share with email

Email a colleague

Share on social media

American Society for Microbiology ("ASM") is committed to maintaining your confidence and trust with respect to the information we collect from you on websites owned and operated by ASM ("ASM Web Sites") and other sources. This Privacy Policy sets forth the information we collect about you, how we use this information and the choices you have about how we use such information.
FIND OUT MORE about the privacy policy