Cross-domain interactions confer stability to benthic biofilms in proglacial streams

Introduction

Biofilms represent the dominant microbial lifestyle in streams and rivers (Battin et al., 2016). These matrix-enclosed microbial communities colonise sediment surfaces, alongside the surface of benthic rocks and plants, and regulate critical ecosystem processes (Battin et al., 2016). Stream biofilm communities are highly diverse, harbouring members from all domains of life, including prokaryotes and various microeukaryotes. This biodiversity enclosed within streams fosters interactions, such as those between bacterial heterotrophs, algae, and fungi. The bacterial and eukaryotic photoautotrophic interactions, the latter of which includes algae, are likely due to metabolic exchanges, whereby algal exudates serve as a source of organic matter for bacterial heterotrophs (Kaplan and Bott, 1989; Wagner et al., 2017). The fungi, meanwhile, may play important roles where fungal metabolites are utilised by distinct bacterial taxa such as Burkholderia (Stopnisek et al., 2016). Furthermore, they are also capable of parasitizing algae, in the case of Chytridiomycetes, to release carbon for bacterial utilisation (Klawonn et al., 2021).

In plants and animals, elevated biodiversity and related biotic interactions can contribute to community stability (McCann, 2000). Yet, studying biotic interactions within complex microbial communities is not a trivial task. For instance, signalling molecules have been suggested as a proxy for identifying interactions among microorganisms (Braga et al., 2016). However, such approaches have remained largely limited to model systems (Vetsigian et al., 2011; Lee and Zhang, 2015), e.g., communities of Streptomyces, Pseudomonas and Vibrio spp. Given the multitude of interacting taxa and the small spatial scales at which interactions occur, the direct observation and quantification of individual microorganisms in stream biofilms is not possible to this date. Instead, patterns of taxa co-occurrence across samples can be used to infer microbial interactions, niches, and key taxa (Faust and Raes, 2012; Berry and Widder, 2014). Despite limitations stemming from the notion that co-occurrence networks are non-empirical and derived from correlations, co-occurrence patterns are often useful for assessing ecological networks and obtaining insights into the organisation of microbial communities (Faust and Raes, 2012; Layeghifard et al., 2017). Co-occurrence networks further allow the exploration of emergent properties, such as the density of interactions, clusters of interacting taxa or the stability of networks against fragmentation. For example, studying bacterial co-occurrences across a dendritic stream network, Widder et al. found evidence for the role of spatial and hydrological processes in shaping co-occurrence network structure and stability (Widder et al., 2014). More recently, Ma et al. highlighted the interconnected patterns across microbiomes in various environments, emphasising the critical impact of microbial interactions on community assembly processes (Ma et al., 2020).

The environment of glacier-fed streams (GFS) is harsh. Low water temperature, coupled with high turbidity, oligotrophy, and snow- and ice-cover over extended periods of the year collectively contribute to this harshness. Unstable stream channels resulting from high turbulence and sediment loads further contribute to challenging conditions, making it difficult for benthic biofilms to establish (Roncoroni et al., 2023). Depending on local topography and geomorphology, GFS may develop into extensive floodplains, with dynamic braided channels exhibiting channel migration on a diel basis (Roncoroni et al., 2023). As distance increases downstream, harshness tends to decrease as channels consolidate and pioneering vegetation stabilises channels (Miller and Lane, 2019). Channel stabilisation increases the habitability of GFS, making these ecosystems more hospitable for benthic biofilms (Roncoroni et al., 2023). Furthermore, towards the edges of the proglacial floodplain, non-glacial tributary streams (TRIB), predominantly fed by groundwater and snowmelt, drain adjacent elevated terrasses that ultimately discharge into the glacier-fed mainstem (Brown et al., 2007). Given the lack of connectivity with the glacier, the TRIB environment is therefore more stable compared to GFS (Roncoroni et al., 2023), with potential consequences for biofilm structure and function (Freimann et al., 2013; Michoud et al., 2023). In fact, despite their close spatial proximity, GFS and TRIB streams host biofilms that differ in their biomass, composition, and diversity, as well as functional potential (Freimann et al., 2013; Brandani and Busi, 2022). It is predicted that as glaciers shrink, GFS will gradually transition into streams fed by groundwater and snowmelt, hence becoming more alike TRIB (Milner et al., 2017; Kohler et al., 2022). However, the consequences of this environmental change for microbial life remain poorly understood.

Brandani et al. recently reported on the complexity of GFS benthic biofilms in the proglacial floodplains of the Swiss Alps (Brandani and Busi, 2022). For instance, alongside bacteria, eukaryotic photoautotrophs such as Ochrophytes (algae) are a major component of these streams. Additionally, fungi such as Chytridiomycetes play key roles in microbial interactions, being involved in 12% of the overall interactions within these streams (Brandani and Busi, 2022). However, the role of taxa within these complex networks in influencing and maintaining community composition have not been resolved. Here, we investigated the properties of cross-domain co-occurrence networks of benthic biofilms in GFS and TRIB with different deglaciation histories in two proglacial floodplains in the Swiss Alps. We focused on putative interactions between bacteria, eukaryotic photoautotrophs, and parasitic fungi. Eukaryotic photoautotrophs are indeed important to stream biofilms as their exudates constitute an energy source for heterotrophic bacteria (e.g. during algal blooms (Kaplan and Bott, 1989), in the presence of light (Battin et al., 2003; Wagner et al., 2017)), while fungi can parasitize algae, thereby affecting microbial carbon flow (Klawonn et al., 2021). Therefore, we posit eukaryotic photoautotrophs and fungi as keystone taxa that take a central role, defined here based on network centrality measures, in maintaining network structure, and further hypothesise that the apparent stability of co-occurrence networks in GFS and TRIB changes along downstream and lateral gradients of deglaciation histories and environmental stability. To evaluate these roles, we assessed the stability of cross-domain co-occurrence networks upon removal of keystone taxa. We also investigated the variance in bacterial community composition that can be explained by specific bacterial and eukaryotic keystone taxa, subsequently contrasting this to the variance that can be explained by environmental differences among sites.

Materials and methods

Sample collection

Benthic biofilms were collected from stream sediments (0 - 5 cm depth) near to the glacier snout and extending to the floodplain’s outlet in various TRIB and GFS stream reaches within the Otemma Glacier (Otemma; 45° 56’ 08.4” N, 7° 24’ 55.1” E) and Val Roseg Glacier (Val Roseg; 46° 24’ 21.1” N, 9° 51’ 55.1” E) floodplains in Switzerland. In each reach, we collected sandy sediments (0.25 - 3.15 mm) with flame-sterilised sieves and spatulas. Sediments samples were transferred to sterile vials, immediately frozen in the dark on dry ice in the field, and stored at -80°C in the laboratory until analysis (all of which were completed within a few months). Samples were collected during early (June/July) and late (August/September) summer (Brandani and Busi, 2022) in 2019. As shown previously, the two sample periods did not show differences in terms of community composition and structure (Brandani and Busi, 2022). For water, i.e. DOC and major ion measurements, samples were refrigerated, while nutrient samples were frozen at -20 C. Samples for the analysis of ions were sterile-filtered and kept at 4C until processing within a few weeks. Samples for the analysis of DOC were GF/F filtered, kept at 4C in the dark and analysed immediately upon return to the laboratory (i.e. within 48h after sampling). Study reaches were categorised into GFS or TRIB depending on their connectivity to glacier runoff based on visual field observations, drone-based imagery, and physicochemical characteristics (Brandani and Busi, 2022). The physicochemical parameters were established previously as part of the study characterising the diversity of the pro- and eukaryotic communities (Brandani and Busi, 2022). Overall, a total of 136 samples (GFS: 50; TRIB: 86) were collected across both floodplains. The sampling sites, indicating pre-2000 (UP) versus post-2000 (DOWN), from GFS and TRIB across Otemma and Val Roseg are depicted in Supplementary Figure 1. These included 68 samples each for the Otemma Glacier and Val Roseg Glacier floodplains, with the exact breakdown of these samples into GFS (Otemma:20; Val Roseg: 30) and TRIB (Otemma:48; Val Roseg: 38) including metadata listed in Supplementary Table 1.

Deglaciation histories

We identified past glacier extents from historic orthophotos and maps using SWISSIMAGE journey through time (Rickenbacher, 2020), and the GLIMS glacier inventory (Raup et al., 2007). These extents were compared with GLAMOS (Linsbauer et al., 2021) frontal variation measurements to verify glacial re-advances. The year of latest glaciation was thus interpolated for each sample site, which provided the longitudinal deglaciation history. We further split the reaches of the floodplain into those which were already deglaciated in the year 2000 (post-2000 or DOWN) and those still glaciated in 2000 (pre-2000 or UP; Supplementary Figure 1A). A lateral gradient (relative to the GFS mainstem) is also given by the TRIB streams that drain the adjacent terraces on the margins of the proglacial floodplains.

Benthic algal biomass

Benthic algal biomass was estimated as chlorophyll α (Supplementary Table 2) using a modified ethanol extraction protocol (Kohler et al., 2020), processed within a few months after collection. For this, sediment (ca. 2 g) samples were treated with 5 ml of 90% EtOH and then placed in a hot water bath (78°C, 10 min), followed by an incubation in the dark (4°C, 24 h). They were thereafter vortexed, centrifuged, and the supernatant read on a plate reader at 436/680 nm (excitation/emission). Chlorophyll α concentrations were inferred from a spinach standard and normalised by the sediment dry mass (DM).

Metabarcoding library preparation and sequencing

A previously established protocol (Busi et al., 2020) utilising phenol-chloroform was used for DNA extraction from benthic sediments (ca. 0.5 g). After initial processing, the DNA samples were diluted to a final concentration of ≤ 2-3 ng/µl. For the 16S rRNA gene metabarcoding analyses, we used the methodology previously described in Fodelianakis et al. (Fodelianakis et al., 2022), where the V3-V4 hypervariable region of the 16S rRNA gene were targeted with the 341F/785R primers. This was done in line with the 16S library preparation Illumina guidelines for the MiSeq system. The eukaryotic 18S rRNA gene metabarcoding library preparation was performed similarly but using the TAReuk454F-TAReukREV3 primers (Stoeck et al., 2010). Based on the MiSeq manufacturer’s protocol, amplicon libraries were prepared where a second PCR was used to add dual indices to the purified amplicon PCR products. This allowed for extensive multiplexing of samples on a single sequencing lane of the MiSeq (Illumina) platform after quantification and normalisation. Samples were subsequently sequenced using a 300-base paired-end protocol in the Lausanne Genomic Technologies Facility (Switzerland).

Metabarcoding analyses

For the 16S and 18S rRNA metabarcoding data analyses, a combination of Trimmomatic v0.36 (Bolger et al., 2014) and QIIME2 v.2020.8 (Bolyen et al., 2019) were used with the latest SILVA database (Quast et al., 2013) v138.1 for taxonomic classification of the gene amplicons, i.e. 16S rRNA and 18S rRNA. From the 16S rRNA amplicon dataset, non-bacterial amplicon sequence variants (ASVs), i.e., archaea, chloroplasts, and mitochondria, were removed from all downstream analyses. The dataset was not rarefied for the analyses owing to the saturation of sequencing curves (Brandani and Busi, 2022) and also based on our previous work where rare taxa play important roles in GFS community structure (Wilhelm et al., 2014). Additionally, the rationale behind discarding the archaeal reads was that the primers used were not designed, and are therefore not optimal, for detecting all lineages of archaea (Bahram et al., 2019). A total of 192 sample libraries were generated for the 16S rRNA sequencing and paired-end sequencing produced a total of 17,200,512 reads, with an average of 89,586 reads per sample. For the 18S rRNA amplicon dataset, a total of 157 amplicon sequence libraries from sediment samples were generated. Meanwhile, singletons and ASVs observed only once were discarded. The paired end 18S rRNA sequencing generated a total of 10,837,518 reads, with an average of 64,127 reads per sample. The 18S ASVs were further clustered into operational taxonomic units (OTUs) based on a 97% identity threshold using the de novo clustering method in vsearch, which has been implemented in QIIME2. This approach was used in order to avoid the overestimation of diversity driven by a high copy number of 18S in eukaryotic cells (Brandani and Busi, 2022). For the analyses, photoautotrophic eukaryotes along with fungi and protists were retained in the 18S rRNA amplicon dataset, given their contributions to community composition as highlighted in our previous work (Battin et al., 2016). Consequently, other non-phototrophic eukaryotes were discarded from the downstream analyses. The 18S rRNA dataset was also not rarefied, and any singletons/OTUs observed in only one sample were removed from downstream analyses, resulting in a dataset of 18S rRNA eukaryotic photoautotrophs and fungi with 429 OTUs. For our downstream analyses (i.e., co-occurrence), it was critical to include paired datasets. This was also a limitation of SpiecEasi (Kurtz et al., 2015), where only paired samples can be included to run pro- and eukaryotic network analyses. To address this, only 136 samples with paired 16S and 18S data were used for downstream analyses.

Co-occurrence networks

To study potential interactions between bacteria, eukaryotic photoautotrophs, and fungi, co-occurrence network analyses were performed with samples meeting specific criteria. These included: 1) the presence of both 16S and 18S sequence data for each sample, and 2) samples had to be categorised the same way across both samplings to ensure replicability (i.e., either designated as GFS or TRIB for both samplings as described by Brandani et al. (Brandani and Busi, 2022). Due to the dynamic nature of proglacial streams, GFS channels tend to meander and migrate, leaving some sites dry, under the influence of TRIB, or infiltrating TRIB streams. Hence, this approach was adopted to avoid potential confounders.

Subsequently, to reduce the noise and overall computational effort, any ASVs found in less than 5% of the samples were discarded from the 16S dataset for the co-occurrence networks. Prior to setting a threshold of 5%, the total number of ASVs retained in the datasets at 1, 5 and 10% thresholds was assessed to reduce noise and computational burden. From a total of 25,307 and 26,912 ASVs in Otemma and Val Roseg, respectively, only 5,268 and 5,216 ASVs were retained at 10%. To avoid exclusivity to only highly-prevalent taxa, we chose the 5% cutoff, where, 17,825 and 16,512 ASVs were respectively retained in Otemma and Val Roseg (Supplementary Table 3).

Co-occurrence networks between 16S and 18S (i.e., bacteria, eukaryotic photoautotrophs and fungi) were constructed using an ensemble of the distance matrices created from SparCC (Friedman and Alm, 2012), Spearman’s correlation (Hauke and Kossowski, 2011), and SpiecEasi (Kurtz et al., 2015) where the networks were constructed using the Meinshausen and Bühlmann (mb) method (Meinshausen and Bühlmann, 2006). Networks were constructed across reaches, for pre-2000 and post-2000segments separately, and across GFS and TRIB for both Otemma and Val Roseg floodplains. Since our analyses are based on amplicon sequence data alone, we focused on the positive interactions across domains to assess potential mutualism within the microbiome. While reports suggest that negative interactions are indicative of co-exclusion mechanisms, especially in human microbiomes (Faust et al., 2012), the paucity of information available, especially in poorly characterised ecosystems may be insufficient to establish via amplicon sequencing data.

To detect communities in the network analyses, we used the Louvain clustering algorithm (Ghosh et al., 2018), removing clusters with less than five nodes. Herein, each community is defined as nodes within the graph with a higher probability of being connected to each other than to the rest of the network. Following this, we calculated network topology measures, including number of nodes and edges, number of clusters, diameter, edge-density, and modularity. These refer to the overall structure of the network such as the number of constituent entities (nodes), their connections (edges), overall length of connections (diameter), and independent cliques (clusters). The adjacency matrix was visualised using the igraph package (Csardi and Nepusz, 2006) in R v4.0.3 (R Core Team, 2021). Centrality measures, degree and betweenness, were also estimated per node, using the igraph v1.3.4 package. The fragmentation (f) of the network was determined as the percentage of the number of disconnected subgraphs over the overall nodes in each network (Widder et al., 2014). Fragmentation was estimated iteratively by the removal of each keystone taxa, i.e., top 10 nodes with both a high degree and a high betweenness in each graph. This information was further used for the subsequent generation of network topologies such as the number of clusters following the initial Louvain clustering of the network.

Community analyses

To explore the role of keystone taxa in structuring biofilm communities, we used constrained ordinations (db-RDA, R function vegan:capscale) using Bray-Curtis distances. In contrast to unconstrained ordination, constrained, or canonical ordination, resolves only the variation that can be explained by the constraints. We first employed a forward selection strategy (vegan:ordistep) to identify sets of non-redundant and significant (p < 0.01) bacterial and eukaryotic keystone taxa that we used as constraints for RDA. We subsequently used the explained variance in the constrained ordination as a measure of how important bacterial and eukaryotic keystone taxa are for explaining bacterial community composition in the two floodplains. Prior to db-RDA, Wisconsin-double standardisation was applied to the ASV counts. The relative abundances of keystone taxa were inspected for multicollinearity and were provided as constraints for stepwise model creation (using 999 permutations). Model significance was evaluated for each RDA axis and explained variance of the constraints was extracted. To contrast variance in bacterial community composition that could be explained by keystone taxa, we performed a similar analysis using environmental parameters. For this, important environmental parameters including pH, water temperature, specific conductivity, dissolved oxygen (DO), turbidity and major ions and nutrients (obtained from Brandani and Busi, 2022) were first standardised (log-transformed) and then supplied to forward selection in db-RDA as described above.

Data analysis

All statistical analyses were performed in R v4.0.3. The ggplot2 (Wickham, 2011) package was used for generating plots in R, while patchwork (https://github.com/thomasp85/patchwork) and Adobe Illustrator were used to arrange the figures as displayed.

Results

Cross-domain interactions underlie stream community structure in proglacial floodplains

In both proglacial floodplains, GFS and TRIB harbour diverse microbial communities which include bacteria, fungi, and phototrophic eukaryotes (Supplementary Figure 1). Previously, we reported on the overall community structure, diversity and abundances of the various microbial components in these streams (Brandani and Busi, 2022), including the role of benthic primary producers in shaping the microbial communities (Michoud et al., 2023). Microalgae such as Chlorophyta, Charophyta, Cryptomonadales, and Ochrophyta form the majority of the phototrophic eukaryotes, which collectively contribute a higher diversity in TRIB compared to GFS (Brandani and Busi, 2022). We previously observed the richness and Pielou’s evenness indices were higher in TRIB compared to GFS. In light of this, we similarly found that bacterial communities were more diverse within the TRIB and dominated by Proteobacteria (Gamma- and Alphaproteobacteria), Bacteroidia, Cyanobacteria, and Verrucomicrobia. Based on covariation of taxa abundances across samples, we built co-occurrence networks. These networks were based on 1,090 nodes including bacteria, eukaryotic photoautotrophs, and fungi, with an average of 61,115 edges (interactions) connecting the nodes. The topological characteristics of the individual networks yielded similar metrics, reflected in density, modularity, assortativity and transitivity (Supplementary Table 1). In all networks, except Otemma post-2000, we identified three dense clusters of co-occurring taxa, one with a majority of eukaryotic photoautotrophs, another comprising mainly bacteria, and an intersecting third cluster composed of microbial eukaryotes including fungi and bacteria (Figure 1).

Figure 1

Figure 1 Network structure of glacier-fed streams and tributary streams. The overall structure of the cross-domain networks from the glacier-fed streams (GFS) and non-glacial tributaries (TRIB) are depicted. (A) GFS from the pre-2000reaches at Otemma, (B) GFS from the pre-2000reaches at Otemma, (C) TRIB from the post-2000 reaches at Otemma, (D) TRIB from the pre-2000reaches at Otemma. From the Val Roseg glacier, the network structures are depicted as follows: (E) GFS from the pre-2000reaches, (F) GFS from the post-2000 reaches, (G) TRIB from the pre-2000reaches, (H) TRIB from the post-2000 reaches. Each node represents a single amplicon sequence variant (ASV), and the lines represent the edges between them, while the colours indicate bacteria, eukaryotic photoautotrophs (indicated as phototrophs) and fungi. The convex hulls indicate clusters identified based on Louvain clustering of the overall network.

Next, we assessed the relative abundance of taxa present in the networks at the family level. Across both floodplains and stream types, we found that Acetobacteraceae were significantly overrepresented in networks constructed from pre-2000compared to post-2000reaches (two-way ANOVA, adj. p < 0.05, Supplementary Figures 2A, B, 3B). On the other hand, Comamonadaceae were significantly overrepresented in post-2000networks (two-way ANOVA, adj. p < 0.05), especially in TRIB (Supplementary Figures 2B and 3B). We also found that Chrysophyceae were overrepresented in pre-2000networks, while Diatomea decreased in pre-2000networks (two-way ANOVA, adj. p < 0.05) (Supplementary Figures 2C, D and 3C, D). Chytridiomycota, parasitic fungi infecting algae (Klawonn et al., 2021), were prevalent in both GFS and TRIB networks, but their abundance did not significantly differ across pre-2000or post-2000sites. However, Zoopagomycota, also parasitic fungi (Spatafora et al., 2016), were considerably enriched in post-2000 reaches across stream types and floodplains (Supplementary Figures 1E-F and 3E, F; adj. p < 0.05, Two-way ANOVA).

Apparent stability of co-occurrence networks

Based on our observations of differential abundance patterns across stream types and deglaciation gradients, we further assessed the contributions of the individual taxa to the overall network. For this, we first identified potential keystone taxa within each network by identifying the top 10 nodes with both a high degree and a high betweenness in each network (Supplementary Figures 4 and 5). For example, taxa classified as Dikarya, Phragmoplastophyta, Chlorophyceae, Cryptomycota, and Diatomea, along with an ASV classified as Burkholderiales, were determined to be keystone taxa in the GFS network at the pre-2000segment of the Otemma Glacier floodplain (Supplementary Figure 4A). Conversely, at the post-2000 segment of the same floodplain, Burkholderiales, Phragmoplastophyta, Xanthophyceae, Chrysophceae, and Dikarya, for instance, were identified as keystone taxa. Similarly, in the pre-2000segment of the Val Roseg Glacier floodplain, Dikarya, Phragmoplastophyta, Gemmatales, Burkholderiales, Cryptomycota, and Diatomea, for instance, were identified as keystone taxa contributing to the network topology (Supplementary Figures 5A, C). Finally, we found various bacteria (e.g., Rhodobacterales, Sphingomonadales) and fungi (e.g., Chytridiomycota) to be keystone taxa in the post-2000 reaches within the Val Roseg floodplain (Supplementary Figures 5B, 5D).

To further understand the role of the keystone taxa in community structure and their effect on apparent network stability, we first assessed network fragmentation upon their removal. For this, the numbers of clusters based on Louvain clustering were determined for each network, following which, a keystone was removed. The fragmentation (f) of the network was assessed before and after iterative removal of the top 10 keystone taxa. Interestingly, we found that in the Otemma Glacier floodplain (Figure 2A), the fragmentation of the networks constructed from the GFS in the post-2000 reaches, increased upon removal of two to three keystone taxa, while the TRIB fragmentation increased upon removal of five keystone taxa. The pre-2000networks, however, appeared more stable, where fragmentation occurred only upon removal of five or eight keystone taxa. In Val Roseg, especially in TRIB (Figure 2B), the overall fragmentation of the microbial network was higher (f_mean=0.48) compared to GFS (f_mean=0.18) upon removal of four or five keystone taxa.

Figure 2

Figure 2 Keystone removal leads to fragmentation of the network. The change in fragmentation (f) for (A) Otemma and (B) Val Roseg are indicated in the line plots, where f was recalculated after each keystone (i.e., bacteria, eukaryotic photoautotroph or phototroph and fungi) was removed from the network. The size of the symbols indicates the relative abundance of the individual keystone taxa within the 16S or 18S data, respectively.

Finally, we unravelled the role of keystone taxa for biofilm community composition. Constrained ordinations revealed that both bacterial as well as eukaryotic phototrophic keystone taxa can explain a substantial fraction of bacterial community dissimilarity at the floodplain scale (Figure 3). Specifically, the relative abundance of bacterial keystone taxa explained 35.0% and 25.4% of variance in bacterial community similarity in Val Roseg and Otemma, respectively. While eukaryotic keystone taxa appeared particularly important for explaining network stability, they played a minor role in explaining bacterial community composition (i.e., 8.5% and 2.4% of explained variance in Val Roseg and Otemma, respectively). This is surprising, particularly in relation to the variance in bacterial community composition that could be explained by environmental conditions, which accounted for a mere 16.5% and 14.5%, respectively. The retained environmental parameters, including streamwater temperature, nutrients and DOC concentration explain differences among TRIB and GFS bacterial communities.

Figure 3

Figure 3 Bacterial keystone taxa well explain bacterial community composition. Constrained ordination of Val Roseg (A–C) and Otemma (D–F) floodplain samples revealed that bacterial keystone taxa (B, E), as identified by their position in co-occurrence networks explained most of the variance in Bray-Curtis distance-based bacterial community composition. This outweighed the role of key environmental parameters (A, D) and of eukaryotic keystone taxa (C, F).

Discussion

Biotic interactions are a salient property of microbial communities, with evidence of cross-domain interactions reported from various ecosystems, including freshwaters (Sun et al., 2022; Wijewardene et al., 2022), oceans (Rowan-Nash et al., 2019), glaciers (Brown and Jumpponen, 2015), and icesheets (Perini et al., 2023). To date, such interactions have not been studied in proglacial stream biofilms. Our findings suggest that biotic interactions, as inferred from co-occurrence patterns, play a pivotal role in influencing the apparent stability of stream biofilm communities along deglaciation and environmental gradients. Although previous reports showed structural and functional differences of the biofilm communities dwelling in different stream types within proglacial floodplains (Freimann et al., 2013; Brandani and Busi, 2022), we found that the overall network topology was similar between both proglacial floodplains, stream types, and deglaciation gradients. This contrasts with our expectation of successional imprints owing to deglaciation on co-occurrence networks. On the one hand, biotic interactions may be established very early on during community succession in the streams that drain recently deglaciated terrain. Indeed, our sampling design covered the successional timescale of the past 20 (pre-2000sites) and 80 (post-2000) years. Furthermore, functional redundancies across clades may also contribute to the apparent similarity of cross-domain interaction networks. Functionally redundant taxa may transiently occupy the same position in interaction networks and therefore result in similar network topologies. However, additional work will be necessary to relate network topology, taxa position and stability with functional characteristics to substantiate this notion.

Cross-domain networks have the potential to reveal key associations between microbial taxa (Williams et al., 2014). We found that biofilms in GFS and TRIB draining recently deglaciated terrain (i.e., pre-2000sites) had relatively more stable networks compared to the post-2000 sites. This finding suggests that bacterial keystone taxa are important for the apparent stability of the cross-domain interaction networks of biofilms dwelling in these nascent stream ecosystems. Furthermore, our results reveal that not all keystone taxa are typically among the most abundant community members, suggesting that low abundance taxa may also play important roles in stabilising microbial networks, corresponding to the notion of keystone species (Han et al., 2022). Our findings agree with observations from recent reports (Crump et al., 2009; Wilhelm et al., 2014; de Cena et al., 2021) highlighting the role of low-abundance taxa in ecosystem function and structure. For example, de Cena et al. recently hypothesised that low-abundance taxa, albeit in the human microbiome, act as keystone taxa, and might often be more metabolically influential within the community (de Cena et al., 2021). Similarly, Crump et al. (Crump et al., 2009) identified microbial keystone taxa, albeit neither rare nor abundant, that are central to ecosystem-level metabolic activity. Additionally, as Ren et al. (Ren et al., 2021) highlighted, there may be several causes including but not limited to environmental variables and/or the metacommunity structure of tributaries in Otemma. It is plausible that both metacommunity structure such as connectivity, directionality of dispersal, and also if these sites are rather discrete communities or part of a gradient/continuum may be contributing to these differences (Leibold and Mikkelson, 2002). It is also likely that the physical heterogeneity of the floodplains, i.e. tributaries in Otemma were larger compared to Val Roseg, may potentially lead to greater heterogeneity and thus increased fragmentation.

Work on multi-trophic food webs (Stouffer and Bascompte, 2011) and agroecosystems (Pocock et al., 2012) has demonstrated the fragility of ecological networks towards removal of key nodes. Our fragmentation analysis substantiates the notion of keystone taxa and their role for the stability of the cross-domain network. Interestingly, we identified several eukaryotes as keystone taxa, underscoring their relevance for biofilm structure and functioning. In GFS in Central Asia, Ren et al. (Ren and Gao, 2019) reported that fungi form integral components of cross-domain interactions networks, forming more clustered networks less susceptible to disturbances. As highlighted previously for stream biofilms (Wagner et al., 2015; Busi et al., 2022), eukaryotic algae serve as sources of organic matter thereby fuelling phototrophic-heterotrophic interactions. Simultaneously, parasitic fungi also foster the release of organic compounds from algae via the ‘fungal shunt’ (Klawonn et al., 2021). The prevalence of parasitic fungi has been noted previously in GFS (Kohler et al., 2022) and other cryospheric ecosystems (Anesio et al., 2017); our analyses further point to the importance of interactions among parasitic fungi and their algal host in proglacial stream biofilms. Along these lines, Mo et al. (Mo et al., 2021) recently suggested that interactions of microeukaryotes between them in the Lena River continental shelf were more stable compared to that of the estuary, potentially explained by variability in salinity. In contrast, Liu and Jiang reported that bacteria-bacteria interactions dominate co-occurrence networks in coastal sea waters of Antarctica (Liu and Jiang, 2020) and related this to competitive abilities of prokaryotes.

Taken together, the roles of bacterial and eukaryotic keystone taxa in ecological networks and their stability may very much be context dependent. We argue that, likely driven by the provisioning of organic matter to heterotrophs, eukaryotic algae and their fungal counterparts play central roles in biofilm interaction networks. However, we quantified the importance of pro- and eukaryotic keystone taxa to overall bacterial community structure and found that the relative abundance of bacterial keystone taxa could explain much of the bacterial community structure. This points towards a hierarchical structuring of interactions among eukaryotic and bacterial members of the biofilm. While eukaryotic primary producers may directly interact with only some bacterial keystone taxa, these bacterial keystone taxa themselves interact, likely via the exchange of secondary metabolites, with a much larger number of bacteria in the biofilm assemblage (Hessler et al., 2023). Such a hierarchical organisation of interactions is likely sensitive to changes in taxa at the base (i.e., the algal primary producers) whereas functional redundancies may dampen the impacts of taxa replacement. This may be particularly relevant in proglacial streams, where reduced light availability due to suspended particles and substrate instability typically inhibit algal growth. The current retreat of glaciers weakens these controls with potential effects on stream microbial communities.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, PRJNA808857 and https://doi.org/10.5281/zenodo.7524289, zenodo-7524289.

Author contributions

SBB: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Visualization, Validation, Writing – original draft. HP: Conceptualization, Formal analysis, Investigation, Methodology, Writing – original draft, Writing – review & editing. JB: Investigation, Writing – original draft, Writing – review & editing. TK: Writing – original draft, Writing – review & editing. SF: Data curation, Writing – review & editing. PP: Methodology, Writing – review & editing. MB: Data curation, Writing – review & editing. GM: Data curation, Writing – review & editing. LE: Data curation, Writing – review & editing. SL: Conceptualization, Funding acquisition, Resources, Writing – original draft, Writing – review & editing. PW: Conceptualization, Funding acquisition, Project administration, Resources, Supervision, Writing – original draft, Writing – review & editing. TB: Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. Funding was provided by the Swiss National Science Foundation grant (CRSII5_180241) to TB, SL and PW.

Acknowledgments

We would like to acknowledge Kevin Casellini, Nicola Deluigi, Matteo Roncoroni, for their help sampling in the field. We would also like to thank Martina Schön for her help with collecting and interpolating glacial extent. The Hunting and Fishing Office of the Canton of the Grisons gave permission to fly the drone on the floodplain of the Tschierva Glacier.

Conflict of interest

The authors 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.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frmbi.2023.1280809/full#supplementary-material

Supplementary Figure 1 | 16S and 18S community profiles. (A) Bird’s eye-view of the Otemma (left) and Val Roseg (right) floodplains depicting the glacier-fed stream mainstem (GFS) and the branching non-glacial tributaries (TRIB). The dashed line indicates the year 2000, where samples were classified as ‘pre-2000 or UPor ‘post-2000 or DOWN’ site above and below, respectively. (B) Family-level profiles of the top 15 bacteria found in the floodplains across reaches and stream types (GFS and TRIB). (C) Relative abundance of the top 15 eukaryotic photoautotrophs.

Supplementary Figure 2 | Taxa contributing to cross-domain interactions in Otemma. Relative abundance of bacteria found in the cross-domain networks of the (A) GFS and (B) TRIB in Otemma. (C) and (D) show the relative abundance of the eukaryotic photoautotrophs in the GFS and TRIB respectively, while (E) and (F) depict the relative abundance of the fungi in Otemma.

Supplementary Figure 3 | Taxa contributing to cross-domain interactions in Val Roseg. Relative abundance of bacteria found in Val Roseg in the cross-domain networks of the (A) GFS and (B) TRIB. Phototroph relative abundances in the (C) GFS and (D) TRIB. (E) and (F) depict the relative abundance of the fungi in GFS and TRIB in Val Roseg.

Supplementary Figure 4 | Keystonetaxa in Otemma. The keystone taxa for the GFS at the (A) pre-2000and (B) post-2000 reaches are highlighted based on their domain of origin. Keystone taxa in the TRIB at the (C) pre-2000and (D) post-2000reaches from the TRIB are simultaneously shown. The x-axis represents the overall betweenness of the individual taxa, whereas the y-axis indicates the degree centrality.

Supplementary Figure 5 | Keystone taxa in Val Roseg. The keystone taxa for the GFS at the (A) pre-2000and (B) post-2000reaches in Val Roseg are highlighted. keystone taxa in the tributaries at the (C) pre-2000and (D) post-2000reaches from the TRIB are depicted in the scatter plots. The x-axis represents the overall betweenness of the individual taxa, whereas the y-axis indicates the degree centrality, i.e., number of connections per node.

Supplementary Table 1 | Metadata and network topology. Glacier metadata, including the glacier from which samples were collected, pre-2000or post-2000 reaches, and type of stream, i.e., Glacier-fed (GFS) or tributaries (TRIB), are indicated along with network topology measures. The dashed line (- - - -) indicates the ‘millennium cut’ based on which samples were classified as ‘pre-2000 or up’ or ‘post-2000 or down’. The solid lines represent the deglaciated history based on the Glacier Extent Database to determine the date since ‘last glaciation’.

Supplementary Table 2 | Chlorophyll-α measurements. Levels of chlorophyll-α measured at the site for each sample are listed along with metadata.

Supplementary Table 3 | Filtering thresholds for network analyses. ASVs were filtered at various thresholds, i.e. 1%, 5% and 10%, to reduce computational burden without losing diversity. Table indicates ASVs retained at chosen thresholds.

Supplementary Table 4 | Pro- and eukaryotic keystone taxa. Abundance information for all ASVs detected in the bacteria (16S) and eukaryote (18S) datasets are provided alongside their indication of Keystone taxa or otherwise.

References

Anesio A. M., Lutz S., Chrismas N. A. M., Benning L. G. (2017). The microbiome of glaciers and ice sheets. NPJ Biofilms Microbiomes 3, 10. doi: 10.1038/s41522-017-0019-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Bahram M., Anslan S., Hildebrand F., Bork P., Tedersoo L. (2019). Newly designed 16S rRNA metabarcoding primers amplify diverse and novel archaeal taxa from the environment. Environ. Microbiol. Rep. 11, 487–494. doi: 10.1111/1758-2229.12684

PubMed Abstract | CrossRef Full Text | Google Scholar

Battin T. J., Besemer K., Bengtsson M. M., Romani A. M., Packmann A. I. (2016). The ecology and biogeochemistry of stream biofilms. Nat. Rev. Microbiol. 14, 251–263. doi: 10.1038/nrmicro.2016.15

PubMed Abstract | CrossRef Full Text | Google Scholar

Battin T. J., Kaplan L. A., Newbold J. D., Cheng X., Hansen C. (2003). Effects of current velocity on the nascent architecture of stream microbial biofilms. Appl. Environ. Microbiol. 69, 5443–5452. doi: 10.1128/AEM.69.9.5443-5452.2003

PubMed Abstract | CrossRef Full Text | Google Scholar

Berry D., Widder S. (2014). Deciphering microbial interactions and detecting keystone species with co-occurrence networks. Front. Microbiol. 5, 219. doi: 10.3389/fmicb.2014.00219

PubMed Abstract | CrossRef Full Text | Google Scholar

Bolger A. M., Lohse M., Usadel B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. doi: 10.1093/bioinformatics/btu170

PubMed Abstract | CrossRef Full Text | Google Scholar

Bolyen E., Rideout J. R., Dillon M. R., Bokulich N. A., Abnet C. C., Al-Ghalith G. A., et al. (2019). Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat. Biotechnol. 37, 852–857. doi: 10.1038/s41587-019-0209-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Braga R. M., Dourado M. N., Araújo W. L. (2016). Microbial interactions: ecology in a molecular perspective. Braz. J. Microbiol. 47 Suppl 1, 86–98. doi: 10.1016/j.bjm.2016.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Brandani P., Busi (2022). Spatial patterns of benthic biofilm diversity among streams draining proglacial floodplains. Front. Microbiol. 13, 948165. doi: 10.3389/fmicb.2022.948165

PubMed Abstract | CrossRef Full Text | Google Scholar

Brown S. P., Jumpponen A. (2015). Phylogenetic diversity analyses reveal disparity between fungal and bacterial communities during microbial primary succession. Soil Biol. Biochem. 89, 52–60. doi: 10.1016/j.soilbio.2015.06.025

CrossRef Full Text | Google Scholar

Brown L. E., Milner A. M., Hannah D. M. (2007). Groundwater influence on alpine stream ecosystems. Freshw. Biol. 52, 878–890. doi: 10.1111/j.1365-2427.2007.01739.x

CrossRef Full Text | Google Scholar

Busi S. B., Bourquin M., Fodelianakis S., Michoud G., Kohler T. J., Peter H., et al. (2022). Genomic and metabolic adaptations of biofilms to ecological windows of opportunity in glacier-fed streams. Nat. Commun. 13, 2168. doi: 10.1038/s41467-022-29914-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Busi S. B., Pramateftaki P., Brandani J., Fodelianakis S., Peter H., Halder R., et al. (2020). Optimised biomolecular extraction for metagenomic analysis of microbial biofilms from high-mountain streams. PeerJ 8, e9973. doi: 10.7717/peerj.9973

PubMed Abstract | CrossRef Full Text | Google Scholar

Crump B. C., Peterson B. J., Raymond P. A., Amon R. M. W., Rinehart A., McClelland J. W., et al. (2009). Circumpolar synchrony in big river bacterioplankton. Proc. Natl. Acad. Sci. U. S. A. 106, 21208–21212. doi: 10.1073/pnas.0906149106

PubMed Abstract | CrossRef Full Text | Google Scholar

Csardi G., Nepusz T. (2006). The igraph software package for complex network research. InterJournal Complex Syst. 1695, 1–9.

Google Scholar

de Cena J. A., Zhang J., Deng D., Damé-Teixeira N., Do T. (2021). Low-abundant microorganisms: the human microbiome’s dark matter, a scoping review. Front. Cell. Infect. Microbiol. 11, 689197. doi: 10.3389/fcimb.2021.689197

PubMed Abstract | CrossRef Full Text | Google Scholar

Faust K., Raes J. (2012). Microbial interactions: from networks to models. Nat. Rev. Microbiol. 10, 538–550. doi: 10.1038/nrmicro2832

PubMed Abstract | CrossRef Full Text | Google Scholar

Faust K., Sathirapongsasuti J. F., Izard J., Segata N., Gevers D., Raes J., et al. (2012). Microbial co-occurrence relationships in the human microbiome. PloS Comput. Biol. 8, e1002606. doi: 10.1371/journal.pcbi.1002606

PubMed Abstract | CrossRef Full Text | Google Scholar

Fodelianakis S., Washburne A. D., Bourquin M., Pramateftaki P., Kohler T. J., Styllas M., et al. (2022). Microdiversity characterizes prevalent phylogenetic clades in the glacier-fed stream microbiome. ISME J. 16, 666–675. doi: 10.1038/s41396-021-01106-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Freimann R., Bürgmann H., Findlay S. E. G., Robinson C. T. (2013). Bacterial structures and ecosystem functions in glaciated floodplains: contemporary states and potential future shifts. ISME J. 7, 2361–2373. doi: 10.1038/ismej.2013.114

PubMed Abstract | CrossRef Full Text | Google Scholar

Friedman J., Alm E. J. (2012). Inferring correlation networks from genomic survey data. PloS Comput. Biol. 8, e1002687. doi: 10.1371/journal.pcbi.1002687

PubMed Abstract | CrossRef Full Text | Google Scholar

Ghosh S., Halappanavar M., Tumeo A., Kalyanaraman A., Lu H., Chavarrià-Miranda D., et al. (2018). “Distributed louvain algorithm for graph community detection,” in 2018 IEEE International Parallel and Distributed Processing Symposium (IPDPS). 885–895.

Google Scholar

Han G., Luong H., Vaishnava S. (2022). Low abundance members of the gut microbiome exhibit high immunogenicity. Gut Microbes 14, 2104086. doi: 10.1080/19490976.2022.2104086

PubMed Abstract | CrossRef Full Text | Google Scholar

Hauke J., Kossowski T. (2011) Comparison of values of Pearson’s and Spearman’s correlation coefficients on the same sets of data. Available at: https://sciendo.com/downloadpdf/journals/quageo/30/2/article-p87.pdf?pdfJsInlineViewToken=1302953392&inlineView=true.

Google Scholar

Hessler T., Huddy R. J., Sachdeva R., Lei S., Harrison S. T. L., Diamond S., et al. (2023). Vitamin interdependencies predicted by metagenomics-informed network analyses and validated in microbial community microcosms. Nat. Commun. 14, 4768. doi: 10.1038/s41467-023-40360-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Kaplan L. A., Bott T. L. (1989). Diel fluctuations in bacterial activity on streambed substrata during vernal algal blooms: Effects of temperature, water chemistry, and habitat. Limnol. Oceanogr. 34, 718–733. doi: 10.4319/lo.1989.34.4.0718

CrossRef Full Text | Google Scholar

Klawonn I., Van den Wyngaert S., Parada A. E., Arandia-Gorostidi N., Whitehouse M. J., Grossart H.-P., et al. (2021). Characterizing the “fungal shunt”: Parasitic fungi on diatoms affect carbon flow and bacterial communities in aquatic microbial food webs. Proc. Natl. Acad. Sci. 118, e2102225118. doi: 10.1073/pnas.2102225118

CrossRef Full Text | Google Scholar

Kohler T. J., Fodelianakis S., Michoud G., Ezzat L., Bourquin M., Peter H., et al. (2022). Glacier shrinkage will accelerate downstream decomposition of organic matter and alters microbiome structure and function. Glob. Change Biol. doi: 10.1111/gcb.16169

CrossRef Full Text | Google Scholar

Kohler T. J., Peter H., Fodelianakis S., Pramateftaki P., Styllas M., Tolosano M., et al. (2020). Patterns and drivers of extracellular enzyme activity in New Zealand glacier-fed streams. Front. Microbiol. 11, 591465. doi: 10.3389/fmicb.2020.591465

PubMed Abstract | CrossRef Full Text | Google Scholar

Kurtz Z. D., Müller C. L., Miraldi E. R., Littman D. R., Blaser M. J., Bonneau R. A. (2015). Sparse and compositionally robust inference of microbial ecological networks. PloS Comput. Biol. 11, e1004226. doi: 10.1371/journal.pcbi.1004226

PubMed Abstract | CrossRef Full Text | Google Scholar

Layeghifard M., Hwang D. M., Guttman D. S. (2017). Disentangling interactions in the microbiome: A network perspective. Trends Microbiol. 25, 217–228. doi: 10.1016/j.tim.2016.11.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee J., Zhang L. (2015). The hierarchy quorum sensing network in Pseudomonas aeruginosa. Protein Cell 6, 26–41. doi: 10.1007/s13238-014-0100-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Leibold M. A., Mikkelson G. M. (2002). Coherence, species turnover, and boundary clumping: elements of meta-community structure. Oikos 97, 237–250. doi: 10.1034/j.1600-0706.2002.970210.x

CrossRef Full Text | Google Scholar

Linsbauer A., Huss M., Hodel E., Bauder A., Fischer M., Weidmann Y., et al. (2021). The new Swiss Glacier Inventory SGI2016: From a topographical to a glaciological dataset. Front. Earth Sci. 9. doi: 10.3389/feart.2021.704189

CrossRef Full Text | Google Scholar

Liu Q., Jiang Y. (2020). Application of microbial network analysis to discriminate environmental heterogeneity in Fildes Peninsula, Antarctica. Mar. pollut. Bull. 156, 111244. doi: 10.1016/j.marpolbul.2020.111244

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma B., Wang Y., Ye S., Liu S., Stirling E., Gilbert J. A., et al. (2020). Earth microbial co-occurrence network reveals interconnection pattern across microbiomes. Microbiome 8, 82. doi: 10.1186/s40168-020-00857-2

PubMed Abstract | CrossRef Full Text | Google Scholar

McCann K. S. (2000). The diversity–stability debate. Nature 405, 228–233. doi: 10.1038/35012234

PubMed Abstract | CrossRef Full Text | Google Scholar

Meinshausen N., Bühlmann P. (2006). Variable selection and high-dimensional graphs with the lasso. Ann. Stat. 34, 1436–1462. doi: 10.1214/009053606000000281

CrossRef Full Text | Google Scholar

Michoud P., Kohler T. J., Peter H., Brandani J., Banu Busi S., Battin T. J. (2023). Unexpected functional diversity of stream biofilms within and across proglacial floodplains despite close spatial proximity. Limnol. Oceanogr. doi: 10.1002/lno.12415

CrossRef Full Text | Google Scholar

Miller H. R., Lane S. N. (2019). Biogeomorphic feedbacks and the ecosystem engineering of recently deglaciated terrain. Prog. Phys. Geography: Earth Environ. 43, 24–45. doi: 10.1177/03091333188165

CrossRef Full Text | Google Scholar

Milner A. M., Khamis K., Battin T. J., Brittain J. E., Barrand N. E., Füreder L., et al. (2017). Glacier shrinkage driving global changes in downstream systems. Proc. Natl. Acad. Sci. U. S. A. 114, 9770–9778. doi: 10.1073/pnas.1619807114

PubMed Abstract | CrossRef Full Text | Google Scholar

Mo Y., Peng F., Gao X., Xiao P., Logares R., Jeppesen E., et al. (2021). Low shifts in salinity determined assembly processes and network stability of microeukaryotic plankton communities in a subtropical urban reservoir. Microbiome 9, 128. doi: 10.1186/s40168-021-01079-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Perini L., Gostinčar C., Likar M., Frisvad J. C., Kostanjšek R., Nicholes M., et al. (2023). Interactions of fungi and algae from the Greenland ice sheet. Microb. Ecol. 86, 282–296. doi: 10.1007/s00248-022-02033-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Pocock M. J. O., Evans D. M., Memmott J. (2012). The robustness and restoration of a network of ecological networks. Science 335, 973–977. doi: 10.1126/science.1214915

PubMed Abstract | CrossRef Full Text | Google Scholar

Quast C., Pruesse E., Yilmaz P., Gerken J., Schweer T., Yarza P., et al. (2013). The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 41, D590–D596. doi: 10.1093/nar/gks1219

PubMed Abstract | CrossRef Full Text | Google Scholar

R Core Team. A language and environment for statistical computing (Vienna, Austria: R Foundation for Statistical Computing). Available at: https://ci.nii.ac.jp/naid/20001692429/.

Google Scholar

Raup B., Racoviteanu A., Khalsa S. J. S., Helm C., Armstrong R., Arnaud Y. (2007). The GLIMS geospatial glacier database: A new tool for studying glacier change. Glob. Planet. Change 56, 101–110. doi: 10.1016/j.gloplacha.2006.07.018

CrossRef Full Text | Google Scholar

Ren Z., Gao H. (2019). Ecological networks reveal contrasting patterns of bacterial and fungal communities in glacier-fed streams in Central Asia. PeerJ 7, e7715. doi: 10.7717/peerj.7715

PubMed Abstract | CrossRef Full Text | Google Scholar

Ren Z., Zhang C., Li X., Ma K., Zhang Z., Feng K., et al. (2021). Bacterial communities present distinct co-occurrence networks in sediment and water of the thermokarst lakes in the yellow river source area. Front. Microbiol. 12, 716732. doi: 10.3389/fmicb.2021.716732

PubMed Abstract | CrossRef Full Text | Google Scholar

Rickenbacher M. (2020). Journeys through time with the Swiss national map series. Proc. 26th Int. Cartographic.

Google Scholar

Roncoroni M., Mancini D., Miesen F., Müller T., Gianini M., Ouvry B., et al. (2023). Decrypting the stream periphyton physical habitat of recently deglaciated floodplains. Sci. Total Environ. 867, 161374. doi: 10.1016/j.scitotenv.2022.161374

PubMed Abstract | CrossRef Full Text | Google Scholar

Rowan-Nash A. D., Korry B. J., Mylonakis E., Belenky P. (2019). Cross-domain and viral interactions in the microbiome. Microbiol. Mol. Biol. Rev. 83. doi: 10.1128/MMBR.00044-18

CrossRef Full Text | Google Scholar

Spatafora J. W., Chang Y., Benny G. L., Lazarus K., Smith M. E., Berbee M. L., et al. (2016). A phylum-level phylogenetic classification of zygomycete fungi based on genome-scale data. Mycologia 108, 1028–1046. doi: 10.3852/16-042

PubMed Abstract | CrossRef Full Text | Google Scholar

Stoeck T., Bass D., Nebel M., Christen R., Jones M. D. M., Breiner H.-W., et al. (2010). Multiple marker parallel tag environmental DNA sequencing reveals a highly complex eukaryotic community in marine anoxic water. Mol. Ecol. 19 Suppl 1, 21–31. doi: 10.1111/j.1365-294X.2009.04480.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Stopnisek N., Zühlke D., Carlier A., Barberán A., Fierer N., Becher D., et al. (2016). Molecular mechanisms underlying the close association between soil Burkholderia and fungi. ISME J. 10, 253–264. doi: 10.1038/ismej.2015.73

PubMed Abstract | CrossRef Full Text | Google Scholar

Stouffer D. B., Bascompte J. (2011). Compartmentalization increases food-web persistence. Proc. Natl. Acad. Sci. U. S. A. 108, 3648–3652. doi: 10.1073/pnas.1014353108

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun P., Wang Y., Huang X., Huang B., Wang L. (2022). Water masses and their associated temperature and cross-domain biotic factors co-shape upwelling microbial communities. Water Res. 215, 118274. doi: 10.1016/j.watres.2022.118274

PubMed Abstract | CrossRef Full Text | Google Scholar

Vetsigian K., Jajoo R., Kishony R. (2011). Structure and evolution of Streptomyces interaction networks in soil and in silico. PloS Biol. 9, e1001184. doi: 10.1371/journal.pbio.1001184

PubMed Abstract | CrossRef Full Text | Google Scholar

Wagner K., Bengtsson M. M., Findlay R. H., Battin T. J., Ulseth A. J. (2017). High light intensity mediates a shift from allochthonous to autochthonous carbon use in phototrophic stream biofilms. J. Geophysical Research: Biogeosciences 122, 1806–1820. doi: 10.1002/2016JG003727

CrossRef Full Text | Google Scholar

Wagner K., Besemer K., Burns N. R., Battin T. J., Bengtsson M. M. (2015). Light availability affects stream biofilm bacterial community composition and function, but not diversity. Environ. Microbiol. 17, 5036–5047. doi: 10.1111/1462-2920.12913

PubMed Abstract | CrossRef Full Text | Google Scholar

Wickham H. (2011). ggplot2: ggplot2. Wiley Interdiscip. Rev. Comput. Stat. 3, 180–185.

Google Scholar

Widder S., Besemer K., Singer G. A., Ceola S., Bertuzzo E., Quince C., et al. (2014). Fluvial network organization imprints on microbial co-occurrence networks. Proc. Natl. Acad. Sci. U. S. A. 111, 12799–12804. doi: 10.1073/pnas.1411723111

PubMed Abstract | CrossRef Full Text | Google Scholar

Wijewardene L., Wu N., Fohrer N., Riis T. (2022). Epiphytic biofilms in freshwater and interactions with macrophytes: Current understanding and future directions. Aquat. Bot. 176, 103467. doi: 10.1016/j.aquabot.2021.103467

CrossRef Full Text | Google Scholar

Wilhelm L., Besemer K., Fasching C., Urich T., Singer G. A., Quince C., et al. (2014). Rare but active taxa contribute to community dynamics of benthic biofilms in glacier-fed streams. Environ. Microbiol. 16, 2514–2524. doi: 10.1111/1462-2920.12392

PubMed Abstract | CrossRef Full Text | Google Scholar

Williams R. J., Howe A., Hofmockel K. S. (2014). Demonstrating microbial co-occurrence pattern analyses within and between ecosystems. Front. Microbiol. 5, 358. doi: 10.3389/fmicb.2014.00358

PubMed Abstract | CrossRef Full Text | Google Scholar