Ecological networks reveal contrasting patterns of bacterial and fungal communities in glacier-fed streams in Central Asia

Study area

The Tian Shan mountains, known as the “water tower of Central Asia,” span a large area of Central Asia from northwestern China to southeastern Kazakhstan and from Kyrgyzstan to Uzbekistan (Fig. 1). In Tian Shan area, glaciers contribute considerably to water resources and play an important role in streamflow regimes (Sorg et al., 2012; Unger-Shayesteh et al., 2013). In China’s portion of the Tian Shan Mountains, there are 7,934 glaciers with a total area of 7,179 km² and a total volume of 707 km³ (Guo et al., 2014). However, glaciers in the Tian Shan Mountains are extremely sensitive to global warming (Kraaijenbrink et al., 2017) and have been shrinking rapidly due to climate warming since the 1970s (Narama et al., 2010; Farinotti et al., 2015). For example, Urumqi Glacier No. 1 (GN1, 43°06′N, 86°49′E) is located in eastern Tian Shan Mountain at the headwater of the Urumqi River (Fig. 1) and retreated from an area of 1.95 km² in 1962 to 1.65 km² in 2010 (Zhang et al., 2014). In 2050, GN1 will likely lose up to 54% of the glacier area and 79% of the ice volume relative to 1980 (Gao et al., 2018).

Field sampling

In June 2016, we investigated two glacier-fed streams in the Tian Shan Mountain. Water samples and benthic biofilm samples were collected from 11 sample sites in total spanning from the elevation of 3,828 to 2,646 m (Fig. 1). The sample sites were chosen along the streams in order to have heterogeneous environments, including different land vegetation and hydrological properties such as glacier contributions to stream flow. However, due to the constraint of accessibility, the sites were not spaced at equal intervals. At each sample site, six to nine submerged rocks were randomly sampled from the stream cross section below 10 cm. A sterilized nylon brush was used to remove the benthic biofilm from each stone in an area of 4.5 cm diameter on the upper surface. The slurry was rinsed with 500 mL sterile water. Approximately 10 mL of the mixed slurry was filtered through a 0.2-μm polycarbonate membrane filter (Whatman, Maidstone, UK) which was immediately frozen in liquid nitrogen in the field. After transported to the lab, the benthic biofilm samples were stored under −80 °C until DNA extraction. In addition, 500 mL water samples were collected for chemical analyses with three replications and stored under 4 °C.

Environmental factors

At each sample site, pH, conductivity (Cond), and elevation were measured in situ using a handheld pH meter (PHI 400 Series; Beckman Coulter, Brea, CA, USA), YSI meter (model 80, Yellow Springs, OH, USA), and GPS unit (Triton 500; Magellan, Santa Clara, CA, USA), respectively. Unfiltered water samples were directly used to measure total nitrogen (TN) and total phosphorus (TP). TN was analyzed by ion chromatography with prior persulfate oxidation (EPA 300.0). TP was analyzed using the ammonium molybdate method with prior oxidation (EPA 365.3). Filtered water samples (filtered with pre-combusted GF/F filters) were used to test nitrate (NO₃⁻), ammonium (NH₄⁺), soluble reactive phosphorus (SRP), and dissolved organic carbon (DOC). NO₃⁻ was analyzed using ion chromatography (EPA 300.0). NH₄⁺ was analyzed using the indophenol colorimetric method (EPA 350.1). SRP was measured according to the ammonium molybdate method (EPA 365.3). DOC was measured using a total organic carbon (TOC) analyzer (TOC-VCPH; Shimadzu Scientific Instruments, Columbia, MD, USA). Water chemistry data was reported in our previous research (Ren, Gao & Elser, 2017).

In glacier-fed streams, both biotic and abiotic environments are tightly linked to the relative contributions of glacier melt and runoff to the stream flow (Milner et al., 2001; Hannah et al., 2007; Kuhn et al., 2011). According to the landscape-based hydrological model proposed by Gao et al. (2016, 2017), we classified the landscape into glaciated and non-glaciated. For each sample site, the proportion of glaciated area (GA) in its sub-catchment was calculated and the proportion of glacier source water (GS) in the total runoff was derived from the model (Gao et al., 2016, 2017). The hydrological distance to glacier terminus (GD) was measured according to the river channel network (Fig. 1). These hydrological parameters were also reported in our previous research (Ren, Gao & Elser, 2017).

In the study area, the vegetation (grassland) status was measured by the normalized difference vegetation index (NDVI) using the Terra moderate resolution imaging spectroradiometer Vegetation Indices (MOD13Q1) Version 6 data downloaded from USGS (https://earthexplorer.usgs.gov/) (Fig. 1). The MOD13Q1 product was generated on June 26, 2017 and has a resolution of 250 m. NDVI is calculated based on the absorption of red light and the reflection of infrared radiation by vegetation (Rouse et al., 1974). The equation is represented as NDVI = (NIR − RED)/(NIR + RED), where NIR is near infrared reflectance and RED is visible red reflectance. It has been demonstrated that NDVI exhibits close relationships with above-ground vegetation biomass and coverage (Carlson & Ripley, 1997; Eastwood et al., 1997; Ren et al., 2019). For each stream site, the average NDVI of its sub-catchment (the upstream area of the stream site) was calculated as the mean NDVI of each pixel in the sub-catchment.

DNA extraction, PCR, and sequencing

Bacterial 16S rRNA gene sequences and fungal 18S rRNA gene sequences were analyzed to determine the bacterial community (BC) and fungal community (FC), respectively. To determine the fungal community, the internal transcribed spacer regions and 18S rRNA genes are commonly used and provide similar results and congruent conclusions (Brown, Rigdon-Huss & Jumpponen, 2014). We used 18S rRNA gene sequencing to detect the fungal community in this study. DNA was extracted using the PowerSoil DNA Isolation Kit (MoBio, Carlsbad, CA, USA) following the manufacturer’s protocol. The V3–V4 regions of the 16S rRNA genes were amplified using the bacterium-specific forward and reverse primers 338F-ACTCCTACGGGAGGCAGCA and 806R-GGACTACHVGGGTWTCTAAT (Invitrogen, Vienna, Austria) (Huws et al., 2007; Masoud et al., 2011; Caporaso et al., 2012). The V4–V5 regions of the 18S rRNA genes were amplified using the fungus-specific forward and reverse primers 817F-TTAGCATGGAATAATRRAATAGGA and 1196R-TCTGGACCTGGTGAGTTTCC-3′ (Invitrogen, Vienna, Austria) (Borneman & Hartin, 2000). The forward primers were barcoded, and the barcodes were designed considering the balanced guanine–cytosine content, minimal homopolymer runs, and no self-complementarity of more than two bases to reduce internal hairpin propensity (Hawkins et al., 2018). PCR reaction systems were prepared using a Premix Taq Kit (Code No. RR902A; Takara, Kusatsu,, Japan) according to the manufacturer’s instructions. The total volume of each PCR reaction was 20 μL, containing 10 μL of 2×EX Premix Taq™ Polymerase, one μL of forward primer, one μL of reverse primer, one μL of NDA extraction, and seven μL of Nuclease-free water. The PCR reactions were conducted with a thermal cycler (ABI 2700, SeqGen, Torrance, CA, USA) with a temperature profile of 1-min hot start at 80 °C, followed by pre-denaturation at 94 °C for 5 min, 30 cycles of amplification (denaturation at 94 °C for 30 s, annealing at 52 °C for 30 s, and extension at 72 °C for 90 s), and a final extension at 72 °C for 10 min. The PCR amplicons were verified in 1.0% agarose with 1× TAE buffer using electrophoresis, purified using the Gel Extraction Kit (Qiagen, Hilden, Germany), and quantified by Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA, USA). One of the fungi sample (N2) was not successfully amplified. The purified and quantified DNA libraries were then pooled together according to their concentrations. The pooled library was sequenced on an Illumina MiSeq (PE300) sequencing platform.

Analyses

Raw sequence data of bacterial 16S rRNA (available at NCBI, PRJNA398147, SRP115356) and fungal 18S rRNA (available at NCBI, PRJNA542974, SRP198430) were processed using QIIME 1.9.0 (Caporaso et al., 2010). The forward and reverse reads were merged. The merged sequences were then assigned to samples based on the barcode. The barcode and primer sequence were cut off to truncate the sequences. The sequences with length >200 bp and mean quality score <20 were discarded. Using UCHIME algorithm (Edgar et al., 2011), the chimeric sequences were detected and removed. Finally, the effective sequences of 16S and 18S rRNA were grouped into OTUs against the Silva 132 database at 99% threshold.

All the analyzed environmental variables, including GD, GA, GS, NDVI, elevation, pH, Cond, TN, NO₃⁻, NH₄⁺, TP, SRP, and DOC, were standardized using the “normalize” method in decostand function in vegan 2.5-3 package (Oksanen et al., 2007). Spearman correlation analyses (cor function in stats v3.6.0 package) were used to assess the pairwise relationships between environmental variables and visualized using corrplot function in corrplot 0.84 package. Mantel tests (mantel function in vegan 2.5-3 package) were used to assess the relationships among spatial dissimilarity (represented by geographic distance), overall environmental distance, hydrological dissimilarity, nutrient dissimilarity, and vegetation dissimilarity. Partial Mantel tests (mantel.partial function in vegan 2.5-3 package) were used to assess those relationships by controlling nutrient dissimilarity. The geographic distance was calculated based on the GPS coordinates of sample sites using distGeo function in geosphere 1.5-10 package. The environmental distance was represented by Euclidean distance based on all analyzed environmental variables. Hydrological dissimilarity was represented by Euclidean distance based on GD, GA, and GS. Nutrient dissimilarity was represented by Euclidean distance based on TN, NO₃⁻, NH₄⁺, TP, and SRP. Vegetation dissimilarity was represented by Euclidean distance based on NDVI. Euclidean distances were calculated using vegdist function in vegan 2.5-3 package (Oksanen et al., 2007). All the above statistical analyses were conducted in R 3.5.1 (R Development Core Team, 2017).

The co-occurrence networks of bacterial and fungal communities were assessed and visualized using Cytoscape (version 3.7.1) (Shannon et al., 2003). In network analyses, the pairwise correlations between OTUs (OTUs with a relative abundance >0.01%) were calculated using Spearman’s correlation based on the relative abundance of OTUs (data was transformed by Hellinger transformation using decostand function in vegan 2.5-3 package). Only strong (R > 0.7 OR R < −0.7) and significant (P < 0.01) correlations were considered in network analysis. ClusterMaker app (Morris et al., 2011) was used to analyze the modular structures of the co-occurrence networks. Modularity values greater than 0.4 suggest that the network has a modular structure (Newman, 2006). The group attributes layout algorithm was used to construct the networks based on modules. The basic topological metrics of networks were calculated, including number of nodes, number of edges, clustering coefficient, characteristic path length, network density, network heterogeneity, and modularity. The taxonomic dissimilarities (beta diversity) of the BC and FC were calculated based on Bray–Curtis distance in terms of the relative abundance of OTUs using R package vegan 2.5-3 (Oksanen et al., 2007). Moreover, the taxonomic dissimilarities of the major modules (modules have more than 15 nodes) were also calculated based on Bray–Curtis distance in terms of the relative abundance of OTUs in the module. BM1, BM2, BM3, and BM4 represent four major modules in the bacterial network. FM1, FM2, FM3, FM4, FM5, and FM6 represent six major modules in the fungal network. Mantel tests were used to assess the correlations between spatial and environmental dissimilarities and taxonomic dissimilarities of bacterial and fungal communities and modules. The relationships among different modules and communities for bacteria and fungi were also assessed using Mantel tests.

Environmental variations and community taxonomic dissimilarities

The environmental variables of the streams varied across the sampling sites and showed differing interrelationships (Fig. 2). GD, GA, GS, NDVI, elevation, pH, and conductivity were closely correlated with each other (Fig. 2A). DOC was negatively correlated with GS, TN, and NO₃⁻, while positively correlated with NDVI and pH (Fig. 2A). Across the sampling sites, hydrological and vegetation dissimilarities had strong linear relationships with overall environmental dissimilarity, while nutrient dissimilarity only had a weak relationship (Fig. 2B). Moreover, hydrological and vegetation dissimilarities had a strong interrelationship with each other but without significant relationships to nutrient dissimilarity (Fig. 2B). Geographic distance only had a significant relationship with nutrient dissimilarity (Fig. 2B). In addition, Mantel tests showed that environmental, hydrological, and vegetation dissimilarities had significant positive relationships to taxonomic dissimilarities of BC but not to FC (Fig. 3). However, geographic distance and nutrient dissimilarity did not have significant relationships with taxonomic dissimilarities of both bacterial and fungal communities (Fig. 3).

Bacteria and fungi co-occurrence patterns and modular structures

In the studied glacier-fed streams, bacteria and fungi formed complex co-occurrence networks (Fig. 4). The bacterial and fungal networks consisted of 904 and 238 nodes and 21,463 and 1,348 edges, respectively (Table 1). The bacterial network had a higher number of nodes and edges as well as a higher network density (Fig. 4; Table 1), indicating the bacterial network was more complex and connected than the fungal network. However, the fungal network exhibited a higher clustering coefficient, characteristic path length, network heterogeneity, and modularity, indicating that the fungal network had a more clustered topology than the bacterial network (Fig. 4; Table 1).

Co-occurrence networks can be compartmentalized into modules within which nodes are closely associated and are expected to share environmental preferences. We found that the OTUs in bacterial and fungal networks were grouped into four and six major modules (modules with more than 15 nodes), respectively (Figs. 4A and 4B). All the modules were formed by various microbial taxa (Figs. 4C and 4D), which shown differently across sample sites (Figs. S1 and S2). Modules had significantly different taxonomic compositions to each other (Tables S1 and S2). Mantel tests showed that the taxonomic dissimilarities of the BC and modules had strong interrelationships (Fig. 5). The taxonomic dissimilarities of FC and modules only had weak interrelationships. Between bacterial and fungal networks, only FM3 had strong relationships with BC and BMs. Mantel tests also revealed that all bacterial modules (BM1, BM2, BM3, and BM4) and one fungal module (FM3) were positively correlated with environmental, hydrological, and vegetation dissimilarities, but not with geographic distance and nutrient dissimilarity (Fig. 3).