Microbial functional genes elucidate environmental drivers of biofilm metabolism in glacier-fed streams

Benthic biofilms in glacier-fed streams harbor diverse microorganisms driving biogeochemical cycles and, consequently, influencing ecosystem-level processes. Benthic biofilms are vulnerable to glacial retreat induced by climate change. To investigate microbial functions of benthic biofilms in glacier-fed streams, we predicted metagenomes from 16s rRNA gene sequence data using PICRUSt and identified functional genes associated with nitrogen and sulfur metabolisms based on KEGG database and explored the relationships between metabolic pathways and abiotic factors in glacier-fed streams in the Tianshan Mountains in Central Asia. Results showed that the distribution of functional genes was mainly associated with glacier area proportion, glacier source proportion, total nitrogen, dissolved organic carbon, and pH. For nitrogen metabolism, the relative abundance of functional genes associated with dissimilatory pathways was higher than those for assimilatory pathways. The relative abundance of functional genes associated with assimilatory sulfate reduction was higher than those involved with the sulfur oxidation system and dissimilatory sulfate reduction. Hydrological factors had more significant correlations with nitrogen metabolism than physicochemical factors and anammox was the most sensitive nitrogen cycling pathway responding to variation of the abiotic environment in these glacial-fed streams. In contrast, sulfur metabolism pathways were not sensitive to variations of abiotic factors in these systems.

help elucidate how natural communities and their functions respond to environmental changes 27,28 . In streams, for example, biofilms are hot spots of microbial activity 29 , contributing substantially to metabolism and biogeochemical cycles through nutrient uptake, transfer of nutrients to higher trophic levels, and remineralization [30][31][32][33] . Unraveling the functional potential of microbial assemblages is critical for gaining insights into their roles in ecosystem processes and their response to environment change 34,35 . For example, microbial metabolism involves a large set of functional genes and biochemical pathways that power biogeochemical cycling in ecosystems. Microbial nitrogen metabolic processes such as nitrification, denitrification, and nitrogen uptake drive nitrogen cycling in streams, controlling nitrogen exports to downstream [36][37][38] . In anaerobic environments, sulfur compounds are major electron carriers for anaerobic microbial metabolism and sulfate reduction is an important process in organic matter degradation [39][40][41] . However, our understanding of microbial functions of benthic biofilm in glacier-fed streams remains extremely limited.
In Central Asia, glacial ecosystems are experiencing unprecedented global warming, hydrological variation, and glacier retreat rapidly [42][43][44][45] . Glacier-fed stream ecosystems in this semi-arid and cold region are particularly vulnerable. Understanding ecosystem functions and predicting the impacts of climate change on glacier-fed stream ecosystems calls for much better knowledge about microbial processes. High throughput molecular technologies and advanced bioinformatics tools make it possible to predict the functional capacity of microbial communities. Following up on our previous study of microbial communities in glacier-fed streams 24 , this research explores the predicted functions of microbial assemblages in benthic biofilms and their relationships to hydrological and physicochemical factors in glacier-fed streams. We predicted metagenomes from 16s rRNA gene sequence data using PICRUSt (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) 46 , identified functional genes associated with nitrogen and sulfur metabolisms based on KEGG (Kyoto Encyclopedia of Genes and Genomes) database 47 , and explored the relationships between metabolic pathways and abiotic factors. We hypothesized that the functional genes of biofilm community in glacier-fed streams are closely associated with community taxonomic compositions and abiotic environments, especially hydrological factors. Although metagenomics reflects potential rather than realized functional capacity, our data offer a window into the poorly known metabolism of benthic biofilms in glacier-fed streams, providing insights into the potential impacts of glacial melting on the ecosystem function.

Methods
Study Area. The Tianshan Mountains are located in Central Asia. Based on the dataset from the Second Glacier Inventory of China 48 , there are 7934 glaciers distributed in the China's Tianshan Mountain with a total area of 7179 km 2 and an ice volume of 707 km 3 , which account for 16.3%, 13.9%, and 15.7% of the total number, total area, and total volume of glaciers in China, respectively 49 . With global warming, these glaciers have been shrinking rapidly [42][43][44][45]  Sampling. In June 2016, we collected water and benthic biofilm samples from 11 sample sites distributed along two glacier-fed streams in the Tianshan Mountains (Fig. 1). The samples from the sites S1-S5 (elevations ranged from 3548 to 2840 m) were collected on June 19, 2016. The samples from the sites N1-N6 (elevations ranged from 3828 to 2646 m) were collected on June 20, 2016. At each sample site, we randomly sampled 6-9 submerged rocks from the river cross section at depths of 10 to 30 cm. The benthic biofilm was removed by rigorously brushing an area of 4.5 cm diameter from the upper surface of each stone using a sterilized nylon brush (changed between samples) and rinsing the slurry with sterile water. From the mixed slurry, approximately 10 ml was filtered through a 0.2-μm membrane filter that was immediately frozen in liquid nitrogen in the field and DNA Extraction, PCR, and Sequencing. Bacterial 16S rRNA genes were analyzed to determine the predicted functional genes and metabolic pathways of the microbial community in biofilms. DNA extraction, PCR, and sequencing were described in our previous study 24 . Physicochemical factors. At each sample site, water temperature (Temp), dissolved oxygen (DO), pH, and conductivity (Cond) were measured in situ using a YSI handheld meter (Model 80, Yellow Springs Instruments, Yellow Springs, Ohio). Elevation was measured using a GPS unit (Triton 500, Magellan, Santa Clara, California). Water samples were collected for nutrients and dissolved organic carbon (DOC) analysis. Total nitrogen (TN) was analyzed by ion chromatography after persulfate oxidation (EPA 300.0) 51 . Nitrate (NO 3 − ) was determined by ion chromatography (EPA 300.0) 51 . Ammonium (NH 4 − ) was analyzed using the indophenol colorimetric method (EPA 350.1) 52 . Total phosphorus (TP) was quantified using the ammonium molybdate method after oxidation (EPA 365.3) 53 . Soluble reactive phosphorus (SRP) was analyzed using the ammonium molybdate method (EPA 365.3) 53 . DOC was analyzed using a Shimadzu TOC Analyzer (TOC-VCPH, Shimadzu Scientific Instruments, Columbia, Maryland). The results of these basic physicochemical factors are shown in our previous study 24 . Hydrological Factors. The glacier distribution data were obtained from the Second Glacier Inventory of China 48 . The river channel network and the catchment boundaries were derived from topographic data, the Digital Elevation Model (DEM) at 90 m resolution (https://earthexplorer.usgs.gov/). Hydrological factors were calculated using the landscape-based hydrological model of Gao et al. (2016Gao et al. ( , 2017 10,17 , including the glaciated area proportion (G A ) for each sample site, runoff proportion generated from glacier source (G S ), and distance to glaciers (G D ). Methods for calculating these basic hydrological factors are shown in our previous study 24 . Analysis. Bacterial 16s rRNA sequence data (available at National Center for Biotechnology Information, SRP115356) were cleaned using the software package QIIME 54 and then clustered to operational taxonomic units (OTUs) with a complete linkage algorithm on a 97% sequence identity level. The metagenomes were predicted from 16s data using Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) 46 . The functional genes associated with nitrogen metabolism and sulfur metabolism were identified based on Kyoto Encyclopedia of Genes and Genomes (KEGG) database 47 . Redundancy analysis (RDA) 55 was conducted to reveal the association of microbial communities in relation to abiotic factors based on functional gene abundances using vegan package 2.4 56 in R 3.3.2 57 . A Monte Carlo test was used to assess the significance of abiotic factors in RDA 55 . Mantel tests were run to assess how functional dissimilarity (Bray-Curtis distance calculated from functional gene abundances) correlated with taxonomic dissimilarity (Bray-Curtis distance calculated from taxon abundances), environmental distance (Euclidean distance), and geographic distance. Correlation analysis was conducted to assess the relationships between metabolic pathways and abiotic factors with the P-value adjusted by Bonferroni correction using psych package 1.7.5 58 in R 3.3.2 57 .

Results and Discussion
Functional Variation. The 16S rRNA data set consisted of 127,146 sequences clustered in 7545 OTUs and 5771 functional genes. The Bray-Curtis distance calculated from functional gene abundances was significantly correlated with operational taxonomic units (OTUs) dissimilarity (Bray-Curtis distance calculated from taxon abundances, Mantel r = 0.614, P < 0.001) and environmental distance (Euclidean distance, Mantel r = 0.567, P < 0.001), but not with geographic distance (Mantel r = 0.287, P = 0.413, Fig. 2), suggesting that the overall changes in functional genes composition were closely associated with community taxonomic composition and abiotic conditions. Redundancy analysis (RDA) was applied to depict the spatial distribution of microbial communities in relation to effects of physicochemical and hydrological factors in terms of the functional gene composition (Fig. 3). The first two axes accounted for 75.6% of the variance (RDA 1: 50.34%; RDA 2: 25.26%). A Monte Carlo test showed that glacier area proportion (G A ), glacier source proportion (G S ), total nitrogen (TN), dissolved organic carbon (DOC), and pH were the abiotic factors most closely associated (p < 0.05) with the distribution of microbial communities (Fig. 3). Canonical correspondence analysis based on taxon abundances showed similar results 24 .
Different biomes clearly harbor distinct microbial assemblages 27,59,60 . It has been demonstrated that metagenomic functional composition is strongly correlated with taxonomic composition in spite of functional redundancy at the genomic level 61,62 . Moreover, recent work also suggests that the metabolic functional potential of microbial communities in ocean and soil is closely related with environmental conditions 63 . Our study also indicates that the functional metabolic potential of glacier-fed streams is strongly influenced by environmental conditions and taxonomic structure. Potential Nitrogen Metabolism. Nitrogen has various chemical species that are cycled by a suite of coupled biogeochemical processes 64 that are catalyzed by various microbe-derived enzymes 65,66 . To deepen our understanding of glacier-fed stream ecosystems, it is crucial to know which N cycling pathways are operating and how they respond to the surrounding abiotic environment. In this study component, we identified and compared the potential nitrogen cycle pathways that might occur in the glacier-fed streams based on the functional genes encoding the enzymes required for specific nitrogen cycling pathways.
Of the 5,571 detected genes, 48 functional genes were associated with nitrogen metabolism ( Figure S1) and the relative abundance of these functional genes ranged from 0.79% to 0.89% of the whole gene reads. Figure S1 summarizes the relative abundance of the genes encoding enzymes that catalyze nitrogen metabolism. The N cycle involves four reduction pathways (assimilatory nitrate reduction, dissimilatory nitrate reduction to ammonia, denitrification, and nitrogen fixation) and two oxidation pathways (nitrification and anammox) (Fig. 4). The relative abundance of genes associated with assimilatory nitrate reduction and dissimilatory nitrate reduction ranged  There are a variety of genes encoding enzymes that catalyze the important transformation reactions of various nitrogen forms ranging from the oxidation states +5 in nitrate to −3 in ammonium. The percentage of associated gene reads were predicted using PICRUSt based on KEGG database. from 0.044% to 0.072% and from 0.066% to 0.151%, respectively (Fig. 5). The assimilatory pathways require energy and start with the reduction of NO 3 − to NO 2 − and then to NH 4 + , which is highly bioavailable and can be readily used by cells for synthesis of amino acids and nucleotides 67,68 . Dissimilatory nitrate reduction to ammonia uses N-compounds to provide energy (ATP) and is another catabolic pathway that can retain the nitrogen in the system in a bioavailable form (NH 4 + ) for further biological processes 69,70 . Generally, assimilatory pathways were much more prevalent than dissimilatory pathways 71 . In our studied glacier-fed streams, however, genes associated with dissimilatory pathways were much more abundant that those involved in assimilatory pathways, suggesting that the microbial assemblages in glacier-fed streams use inorganic nitrogen more as an energy source more as an N source for biosynthesis. Nitrogen fixation is an energetically demanding process and only a few taxonomic groups genera can carry out this process 72,73 . As a competing process to nitrate reduction pathways, denitrification is the main biological process for removal of N from freshwater systems 70,74 . The relative abundance of genes associated with denitrification and nitrogen fixation ranged from 0.038% to 0.123% and from 0.012% to 0.074%, respectively (Fig. 5). Nitrification is an essential process in nitrogen cycle performed by nitrifiers converting ammonium to nitrate. In our study, functional genes encoding enzymes that catalyze nitrification had a relative abundance ranges from 0.01% to 0.05% (Fig. 5). Anammox plays a significant role in the nitrogen cycle, converting nitrite and ammonia to dinitrogen 75 , and is another microbial process that releases fixed nitrogen from the environment as dinitrogen [76][77][78][79] . In our study, functional genes encoding enzymes that catalyze anammox had only a very small relative abundance (Fig. 5).
In ecosystems, nitrogen metabolic pathways are particularly susceptible to environmental redox fluctuations because of the large difference between the oxidation state of nitrate (+5) in nitrate and that of ammonia (−3) 80,81 . Moreover, some other environmental factors also strongly impact nitrogen cycling. For example, microorganisms can be stimulated by high concentration of NO 3 − to increase populations containing the functional genes encoding the enzymes required for nitrate reduction. Low N:P ratios and low NH 4 + can stimulate N 2 fixing bacteria [82][83][84] . It has also been demonstrated that anammox is influenced by physicochemical properties, such as nitrate concentration, C/N ratio, ammonia concentration, and pH 85 . Nitrification can be affected by temperature, salinity, light, organic matter concentration, substrate concentrations, pH, and oxygen concentration 86 . In our study, we explored the relationships between nitrogen cycle pathways and abiotic environment including physicochemical and hydrological factors (Table S1). The results showed that anammox was negatively correlated with G A and G S . Denitrification, dissimilatory nitrate reduction, and nitrification were negatively correlated with G D . In addition, anammox was positively correlated with pH, temperature, and DOC but negatively correlated with TN and NO 3 . Assimilatory nitrogen reduction was positively correlated with TN and NO 3 . Nitrification was negatively correlated with pH. Nitrogen fixation was positively correlated with NH 4 . These results suggest that hydrological factors had more significant influences on nitrogen metabolism than physicochemical factors. Among the nitrogen Figure 5. Relative abundance of major categories of functional genes encoding the enzymes that catalyze nitrogen cycling pathways (dissimilatory nitrate reduction, assimilatory nitrate reduction, denitrification, nitrification, nitrogen fixation, and anammox) based on the KEGG database. Figure 6. Relative abundance of major categories of functional genes encoding the enzymes that catalyze sulfur metabolism, including assimilatory sulfate reduction, dissimilatory sulfate reduction, and the SOX system.
Scientific REpoRts | 7: 12668 | DOI:10.1038/s41598-017-13086-9 cycling pathways, anammox was the most sensitive to variation of abiotic environment in the glacial-fed streams we studied (Table S1). Since global warming leads to significant changes of abiotic environment, especially in the hydrology in glacier-fed streams in Central Asia, these results may allow us to predict changes in nitrogen metabolism under climate change.
Potential Sulfur Metabolism. Overall, there were 43 functional genes associated with sulfur metabolism in our data set ( Figure S2). The relative abundance of these functional genes ranged from 0.848% to 1.009% of the total gene reads. Figure S2 summaries the relative abundance of the genes encoding sulfur metabolism based on the KEGG database. Assimilatory sulfate reduction occurs in many anoxygenic phototrophic bacteria and leads to the biosynthesis of sulfur-containing compounds 87,88 . Dissimilatory sulfate reduction is a pathway in which sulfate or sulfur serves as the terminal electron acceptor of the respiratory chain, producing inorganic sulfide. In our study, the relative abundance of genes associated with assimilatory sulfate reduction ranged from 0.247% to 0.313%, followed by genes involved in the sulfur oxidation (SOX) system and in dissimilatory sulfate reduction (Fig. 6), suggesting that microbial assemblages in glacier-fed streams regard inorganic sulfur as an S source in biosynthesis more than an electron acceptor. The SOX system, a sulfur oxidation pathway, is a multienzyme pathway found in both photosynthetic and lithoautotrophic sulfur-oxidizing bacteria 39,89,90 . Sulfur oxidation is a critical step in the biogeochemical sulfur cycle 91,92 . Our detection of genes encoding the enzymes catalyzing sulfur oxidation suggests that microbial activity is generating acidity that may also affect mineral formation and stability 88 . Previous studies have suggested that the availability of electron donors 93 , light 94 , and aerobic conditions 95 are important environmental factors affecting sulfur metabolism. In our study, only assimilatory sulfate reduction was positively correlated with G D (Table S1), suggesting that sulfur metabolism pathways were not sensitive to variations of hydrological and physicochemical factors in glacier-fed streams.

Conclusions
In glacier-fed streams, microbial assemblages in benthic biofilms change significantly with variation in the abiotic environment and thus appear vulnerable to glacier retreat caused by climate change. As most glaciers on Earth are receding rapidly, the habitat template (especially hydrological factors) of glacier-fed streams will be increasingly disrupted. Community taxonomic composition and the abiotic environment regulate microbe-mediated functions. In particularly, hydrological factors have significant impacts on the distribution of functional genes and nitrogen cycling pathways. Thus, large changes in biogeochemical processes and associated ecosystem functions likely will occur as glacier-fed ecosystems are transformed due to loss of their association with the cryosphere.