TABLE 1 | Estimation of Grand Calumet River as a potential source† of bacterial populations (mean percentages ± standard error of mean) on three sampling dates for (A) Lake communities as sink‡ , or (B) Specific nearshore communities as sink§ .

†Source contribution estimated using SourceTracker (Knights et al., 2011) plugin in QIIME 1.9.1 (Caporaso et al., 2010) using default conditions. ‡Sinks are river mouth, nearshore, and offshore samples a (n = 3, 9, 6 per sampling date, respectively), Krustal Wallis (P = 0.00003). §Sinks are three nearshore beaches: Jeorse Park (JP), Whihala (WH), and 63rd Street (63) in Chicago, Krustal Wallis (P = 0.004). Different superscripts letters in the tables indicate significant difference of pairwise comparisons by date using Dunn's post hoc test with Bonferroni correction.

the phylum Proteobacteria whereas in the JP samples they were mainly in the Bacteroidetes and a mixture of Proteobacteria and Actinobacteria in the site 63rd and WH samples. Bacteroidetes are a common phylum in fecal samples but since the taxa that differed are not classified to the genus level it is not possible to speculate if they are feces associated. However, at LDA 2.0 unclassified Enterobacteriaceae are significant in the JP sample (**Supplementary Table S4**). This group associated with fecal contamination. Of the nearshore site JP had the highest frequency of elevated E. coli MPN values (**Supplementary Table S5**). This indicates that targeted metagenomic analysis can provide additional data of fecal contamination.

### Microbial Community and E. coli

Although the incidences of elevated E. coli MPN values were limited, correlation analysis indicated that there were positive correlations with other biomarkers and bacterial genera commonly reported from fecal samples. There were

few samples that had elevated E. coli (>235 CFU/MPN per 100 ml) that would lead to closures of recreational waters (**Supplementary Table S5**; U.S. EPA, 2012). Spearman Rho analysis of bacterial taxa relative abundances and E. coli MPN (sorted into bins) indicated there were 52 genera with significant positive correlations and 1 negative correlation (**Supplementary Table S6**). Notable is the positive correlation with several genera in the phylum Firmicutes that are associated with fecal contamination, Enterococcus, Blautia, Faecalibacterium, Clostridium, Collinsella, and Selenomonas as well as Bacteroides from the phylum Bacteroidetes.

Enterococcus (Wheeler et al., 2002) and Bacteroides (Bernhard and Field, 2000) are often used as fecal indicator species, and more recently, Blautia (Koskey et al., 2014; Eren et al., 2015) and Dialister (Jeong et al., 2011) have also been proposed as indicator species. A species within Clostridium (C. perfringens) has been suggested as a reliable indicator of water quality in tropical areas (Fujioka and Shizumura, 1985; Fujioka et al., 1997), where traditional indicators such as E. coli, and enterococci, are commonly found in the environment (Byappanahalli and Ishii, 2011). Collinsella, Blautia, and Faecalibacterium, are examples of commonly found members of gut microbiomes that correlated with E. coli MPNs. With the availability of high throughput sequencing technology, others have suggested that community analysis may be an additional means to assess water quality and can be applied to microbial source tracking (Newton et al., 2013; Henry et al., 2016).

The influence of physical conditions on overall microbial community was also examined (**Figure 6**) using CCA. Factors corresponding to the differences in the bacterial communities were water type (river water being most distant from lake) along the first axis explaining 64.3% of the variation; sampling time, water temperatures and dissolved oxygen along the second axis (14.3%) (**Figure 6A**); and E. coli MPN and turbidity along the third and fourth axes (**Figure 6B**) (100 permutations, p < 0.01). The influence of water type describing the microbial community is like the finding of beta diversity analysis due to differences of the river community to the other water types (**Figure 3A**). A secondary factor was the difference in communities with sampling time (**Figure 3B**) that was shown by CCA to be influenced by increasing dissolved oxygen and decreasing water temperatures (**Figures 6A,B**). The E. coli MPN corresponded to the three river samples collected on the third sampling date after rainfall (**Figure 6B**). This suggests that both shoreline sources (e.g., gulls, runoff, nearshore sand-water interactions) and large-scale processes (e.g., waves, currents, lake turnover), as well as time of year, are likely to influence changes in community along the shoreline. The relative impact of these factors could have inter-annual variation. Integration of physical modeling and multiple years of data could help resolve some of these interactions.

### Implications for Mitigation Strategies

After the recognition that aquatic systems have become degraded, it is essential to develop strategies to restore compromised ecosystems through appropriate remedial actions and in many instances, those actions tend to be site-specific. For instance, at the study locations, nonpoint sources of microbial contamination by shoreline birds (gulls), has been previously identified as a major contributing factor using microbial source tracking and gull deterrence activities have significantly improved shoreline water quality in recent years (Nevers et al., 2018). Similarly, the periodic presence of the human marker Bacteroides HF183 (Byappanahalli et al., 2015) has been a cause of concern for the potential significant human health threat; previously, the Grand Calumet River was the presumptive source of any human contamination. Results presented here, however, indicate that the river, even following a rain event, was likely not impacting shoreline bacterial communities. There may, therefore,

be additional sources of human contamination currently unknown impacting these shoreline locations. At Jeorse Park and 63rd, the shoreline configuration has also been identified as contributing to the persistence of bacterial contamination due to the tendency toward accumulation and the decrease in circulation (Ge et al., 2012; Byappanahalli et al., 2015). Because the river is not the major exogenous source of microbes into the nearshore lake sites, efforts to curtail river flow toward recreational beach areas would not be a sufficient means to decrease beach closures. Efforts toward identifying and mitigating shoreline or nonpoint sources, such as decreasing gull presence (Nevers et al., 2018) and reducing other contributions (e.g., shoreline algae Whitman et al., 2014b), would likely have a greater impact on eliminating beach closures. By incorporating emerging technologies, such as the microbiome, into water quality monitoring programs will be helpful to study long-term changes and underlaying factors that influence these changes which are difficult to elucidate with traditional monitoring programs.

To summarize using a targeted 16S rRNA gene sequencing approach to analyze aquatic microbial communities demonstrated it is possible (1) to identify sources of contamination using a consortium of microbes as an index of pollutant sources rather than using one or two traditional indicators (e.g., E. coli, enterococci) often used in monitoring programs. Using this approach, we were able to show that the Grand Calumet River had minimal influence on shoreline water quality at the study sites, indicating that the sources contributing to high E. coli levels (e.g., at JP, WH) were more likely internal (e.g., shoreline birds, sand, Cladophora) with an intermittent contamination from GCR. (2) The data gathered from this research will be useful to management agencies, such as U.S. EPA and Indiana Department of Environmental Management, for addressing water quality restoration efforts currently under implementation at the study locations. An example will be delisting shorelines as impaired for beneficial use such as recreation. This can begin the process to modify acceptable approaches by government agencies to identify microbially contaminated locations to implement effective mitigation.

### AUTHOR CONTRIBUTIONS

All authors wrote and approved the final version of the manuscript. MN and MB conceived the project. MB, MN,

### REFERENCES

and CN performed the sample designing, contributed to methodology, and wrote the manuscript. CN analyzed the data.

### FUNDING

Funding was provided to the USGS through the Great Lakes Restoration Initiative and USGS agreement G15AC00427, and to Purdue University from the USDA National Institute of Food and Agriculture Hatch grant IND010825.

### ACKNOWLEDGMENTS

We thank P. Ryan Jackson, Kasia Kelly, Dawn Shively, Ashley Spoljaric, and Paul Buszka, USGS and Stephanie Grebinoski, Purdue University for the collection of samples and for help during sample processing. The Purdue Genomics Core Facility at Purdue University performed sequencing and the Purdue Bioinformatics Core kindly provided servers to run the bioinformatics analyses.

### SUPPLEMENTARY MATERIAL

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

FIGURE S1 | Rarefaction curves of alpha diversity measures (A) Shannon Diversity and (B) PD whole tree (Faith's Phylogenetic Diversity). Each line represents a single sample.

TABLE S1 | Relative proportion of phyla represented in Lake Michigan water sample.

TABLE S2 | LEfSe results of a comparison of different Lake Michigan water sources, taxa at all levels with LDA > 2.0 are listed.

TABLE S3 | LEfSe results of a comparison of Lake Michigan sampling sites, taxa at all levels with LDA > 2.0 are listed GCR.

TABLE S4 | LEfSe results of a comparison of only the nearshore Lake Michigan sampling sites, taxa at all levels with LDA > 2.0 are listed.

TABLE S5 | Environmental conditions and E. coli most probable numbers (MPN) in triplicate samples from each sampling location on three collection dates.

TABLE S6 | Significant Spearman rho correlations E. coli most probable numbers (MPN) binned and relative abundances of bacteria at genus level.

communities associated with decaying Cladophora in Great Lakes. Sci. Total Environ. 574, 872–880. doi: 10.1016/j.scitotenv.2016.09.107

**Conflict of Interest Statement:** 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.

Copyright © 2019 Nakatsu, Byappanahalli and Nevers. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# The DNRA-Denitrification Dichotomy Differentiates Nitrogen Transformation Pathways in Mountain Lake Benthic Habitats

Carlos Palacin-Lizarbe1 \*, Lluís Camarero2† , Sara Hallin3‡ , Christopher M. Jones3§ , Joan Cáliz2||, Emilio O. Casamayor2¶ and Jordi Catalan1,4#

1 Centro de Investigación Ecológica y Aplicaciones Forestales, Cerdanyola del Vallès, Spain, 2 Center for Advanced Studies of Blanes, (CEAB–CSIC), Girona, Spain, 3 Department of Forest Mycology and Plant Pathology, Swedish University of Agricultural Sciences, Uppsala, Sweden, 4 Consejo Superior de Investigaciones Científicas, Cerdanyola del Vallès, Spain

### Edited by:

Haihan Zhang, Xi'an University of Architecture and Technology, China

### Reviewed by:

Marja Tiirola, University of Jyväskylä, Finland Katharina Kujala, University of Oulu, Finland

### \*Correspondence:

Carlos Palacin-Lizarbe<EMAIL_ADDRESS>orcid.org/0000-0002-1572-6053 †orcid.org/0000-0003-4271-8988 ‡orcid.org/0000-0002-9069-9024 §orcid.org/0000-0002-2723-6019 ||orcid.org/0000-0001-5714-7431 ¶orcid.org/0000-0001-7074-3318 #orcid.org/0000-0002-2934-4013

### Specialty section:

This article was submitted to Aquatic Microbiology, a section of the journal Frontiers in Microbiology

Received: 15 January 2019 Accepted: 16 May 2019 Published: 04 June 2019

### Citation:

Palacin-Lizarbe C, Camarero L, Hallin S, Jones CM, Cáliz J, Casamayor EO and Catalan J (2019) The DNRA-Denitrification Dichotomy Differentiates Nitrogen Transformation Pathways in Mountain Lake Benthic Habitats. Front. Microbiol. 10:1229. doi: 10.3389/fmicb.2019.01229 Effects of nitrogen (N) deposition on microbially-driven processes in oligotrophic freshwater ecosystems are poorly understood. We quantified guilds in the main N-transformation pathways in benthic habitats of 11 mountain lakes along a dissolved inorganic nitrogen gradient. The genes involved in denitrification (nirS, nirK, nosZ), nitrification (archaeal and bacterial amoA), dissimilatory nitrate reduction to ammonium (DNRA, nrfA) and anaerobic ammonium oxidation (anammox, hdh) were quantified, and the bacterial 16S rRNA gene was sequenced. The dominant pathways and associated bacterial communities defined four main N-transforming clusters that differed across habitat types. DNRA dominated in the sediments, except in the upper layers of more productive lakes where nirS denitrifiers prevailed with potential N2O release. Loss as N2 was more likely in lithic biofilms, as indicated by the higher hdh and nosZ abundances. Archaeal ammonia oxidisers predominated in the isoetid rhizosphere and rocky littoral sediments, suggesting nitrifying hotspots. Overall, we observed a change in potential for reactive N recycling via DNRA to N losses via denitrification as lake productivity increases in oligotrophic mountain lakes. Thus, N deposition results in a shift in genetic potential from an internal N accumulation to an atmospheric release in the respective lake systems, with increased risk for N2O emissions from productive lakes.

Keywords: denitrification, DNRA, lithic biofilms, mountain lake, nitrogen deposition, remote ecosystems, sediment, 16S

### INTRODUCTION

According to the planetary boundaries framework (Rockström et al., 2009), anthropogenic alteration of the nitrogen (N) cycle is one of the major challenges facing the Earth system. Human activities have at least doubled the levels of reactive N (Nr) available in the biosphere (Erisman et al., 2011), resulting in deposition of Nr in or near heavily populated areas as well as remote ecosystems (Catalan et al., 2013). In the context of global change, remote ecosystems — defined here as being affected by atmospheric processes rather than direct human action in catchment areas — can be particularly informative about potential large-scale changes in the Earth system

(Catalan et al., 2013). Alpine lakes of the Northern hemisphere and subarctic regions are examples of remote ecosystems that have been exposed to increased Nr deposition during the last decades (Holtgrieve et al., 2011; Camarero, 2017), triggering a nutrient imbalance in these freshwater systems which are otherwise known to have low nutrient availability (Catalan et al., 2006). While alpine and subarctic lakes are often considered important sensors of global change (Smol, 2012), there is minimal understanding of how increased Nr availability affects microbially-driven N-cycle pathways in these ecosystems (McCrackin and Elser, 2010; Palacin-Lizarbe et al., 2018).

The N cycle is best described as a modular and complex network of biological N-transformation reactions carried out by metabolically versatile communities of microorganisms (Graf et al., 2014; Kuypers et al., 2018), whose overall composition largely determines whether Nr is lost, via denitrification or anammox, or retained in the system via dissimilatory nitrate reduction to ammonium (DNRA). Within lakes, benthic habitats are known as hotspots of N cycling due to steep redox gradients in the sediments and biofilms (Melton et al., 2014). Furthermore, the presence and composition of macrophytes also influence the biogeochemistry of the sediment (Gacia et al., 2009). In particular, isoetid species oxygenate the sediment and may promote coupled nitrification-denitrification (Vila-Costa et al., 2016). However, the effect of increased N deposition on the N-cycling microbial communities, and the factors controlling their distribution are poorly understood in mountain lakes.

Our study aims to investigate how the distribution of microbial communities in general and those that drive different N-transformation pathways changes across a range of different benthic habitats in mountain lakes that have been affected by enhanced N deposition in the absence of significant acidification (Camarero and Catalan, 1998). We hypothesise that benthic habitat type and lake productivity together determines the fate of deposited N and that increased productivity will promote pathways resulting in Nr loss. Lakes at lower altitudes tend to be more productive, particularly if they are small since the productive period is longer (Catalan et al., 2009) and phosphorus loading to the lake increases as the catchment is more vegetated (Kopàcek et al., 2011 ˇ ). In the Pyrenees, more than 70% of the lakes are considered ultraoligotrophic based on total phosphorus (TP; <150 nM), whereas 22 and 6% are oligotrophic and mesotrophic, respectively (Catalan et al., 2006). In general, more productive oligotrophic mountain lakes exhibit low dissolved inorganic nitrogen (DIN) concentrations due to higher consumption of excess N from atmospheric loading by primary producers (Camarero and Catalan, 2012). We therefore selected lakes to establish a DIN gradient and sampled lithic biofilms, sediments with elodeid, isoetid and helophyte macrophytes, and littoral and deep non-vegetated sediments (**Figure 1**). We then characterised the N-functional pathways by quantifying the abundances of key N-functional genes involved in denitrification, nitrification, DNRA and anammox pathways (**Table 1**). We also determined the bacterial community composition in the benthic habitats and linked these to the functional guilds using a multivariate approach combined with indicator species analyses. The environment was characterised by including proximal (benthic) and more distal (lake) descriptors to capture potential drivers acting at different spatial scales (Wallenstein et al., 2006; Battin et al., 2016).

# MATERIALS AND METHODS

# Sampling Location and Habitat Description

The lakes are located in the central region of the Pyrenees mountain range within the Aigüestortes i Estany de Sant Maurici National Park (**Table 2** and **Figure 1**). All lakes are dimictic and ultra-oligotrophic (TP < 150 nM) except for Bassa de les Granotes, which is classified as oligotrophic (150 < TP < 300 nM; Catalan et al., 1993) with a circumneutral pH (∼7; Vila-Costa et al., 2014). All main benthic habitats in the lakes were considered (**Figure 1**), although certain habitats were present in only a few lakes (**Table 2** and **Supplementary Table S1**). Plan Lake is particularly rich in macrophytes, including isoetids (Isoetes setacea, I. palustris, and Subularia aquatica), elodeids (Myriophyllum alterniflorum, Potamogeton alpinus, and P. berchtoldii) and the helophyte Carex rostrata (Gacia et al., 1994). Sampling was carried out during the ice-free period (June–November) of 2013 and 2014, with a total of 30 sites and 226 samples analysed.

# Water, Lithic Biofilm, and Sediment Characterisation

The overlying water, sediments and lithic biofilms were characterised using physical, chemical and biological variables (**Supplementary Table S1**). The temperature of the overlying water was measured at the time of sampling. For chemical analyses, water samples were filtered through a pre-combusted (4 h at 450◦C) GF/F glass fibre philtre. Nitrate and sulphate were determined by capillary electrophoresis using a Quanta 4000 (Waters) instrument. Ammonium and nitrite were determined by colourimetric methods in a segmented-flow autoanalyser (AA3HR, Seal), using the Berthelot reaction for ammonium (Bran+Luebbe method G-171-96) and the Griess reaction for nitrite (Bran+Luebbe method G-173-96). Dissolved organic carbon (DOC) was measured by catalytic combustion to CO2 and detection by IR spectroscopy in a TOC5000 (Shimadzu) analyser.

Lithic biofilms were sampled collecting several cobbles (ø ∼10 cm) from different sites of the lake. Cobbles were scraped entirely (upper and lower sides) with clean metal brushes and washed with deionized water and pooling together the collected material. Biofilm subsamples were collected on 0.2-µm pore polycarbonate membranes for DNA analysis, and triplicate volumes were filtered through a pre-combusted and pre-weighted GF/F glass fibre philtre for chemical and physical analyses. Sediment cores (ø 6.35 cm) were collected with a gravity corer (Glew, 1991) around the deepest point of each lake or manually by scuba diving for the littoral sediments. The cores were sliced in three sections (0–0.5, 0.5–2, and 2–4 cm) to capture the oxic and the nitrate reduction zones (Melton et al., 2014).

habitats: sediments near the deepest point of the lake (non-vegetated) (N,O); littoral sediments from beds of isoetids [Isoetes lacustris (T)] and elodeids [Myriophyllum alterniflorum (U), Potamogeton alpinus (V)] macrophytes, helophyte (Carex rostrata) belts (H–L) and rocky areas (B–E); and lithic biofilms from littoral cobbles (P–R).

For total carbon (C) and N and isotopic composition, ca. 5 mg of the freeze-dried sample was placed with a catalyst (Va2O5) in tin capsules, and the analyses were performed by the

University of California Davis Stable Isotope Facility. Organic matter (OM) content was determined using the loss on ignition (LOI) procedure (Heiri et al., 2001). The median grain size of the sediment was determined by laser diffraction (Mastersizer, 2000, Malvern Instruments Ltd, United Kingdom), using freezedried sediment rehydrated in distilled water and introduced into the sample dispersion unit (Hydro 2000 G, Malvern Instruments Ltd, United Kingdom) after adding hexametaphosphate and sonicating to avoid aggregates. Laser obscuration was between 10–20% and the measuring range between 0.02 and 2000 µm.

### DNA Extraction and Quantitative PCR of 16S rRNA and N Cycle Genes

DNA was extracted from 0.33 ± 0.06 g of sediment or lithic biofilm using the FastDNA R Spin Kit for Soil (MP Biomedical) following the manufacturer's instructions. The extracted DNA was quantified using the Qubit R fluorometer (Thermo Fisher Scientific Inc.).

Quantitative real-time PCR (qPCR) was used to quantify functional genes encoding enzyme involved in N-cycle pathways

(**Table 1** and **Supplementary Table S2**), as well as the bacterial 16S rRNA gene. All qPCR reactions were performed in duplicate in a total reaction volume of 20 µL using DyNAmo Flash SYBR Green qPCR kit (Thermo Fisher Scientific Inc.), 0.1% Bovine Serum Albumin, 0.5–1.0 µM of each primer and 15 ng DNA on the Biorad CFX Connect Real-Time System. Primers (Rotthauwe et al., 1997; Hallin and Lindgren, 1999; Lopez-Gutierrez et al., 2004; Mohan et al., 2004; Throbäck et al., 2004; Henry et al., 2006; Schmid et al., 2008; Tourna et al., 2008; Jones et al., 2013; Welsh et al., 2014), amplification protocols and resulting efficiencies for each assay are listed in **Supplementary Table S3**. Potential inhibition of the PCR reactions was checked by amplifying a known amount of the pGEM-T plasmid (Promega) with the plasmid-specific T7 and SP6 primers added to the DNA extracts and non-template controls. No inhibition of the amplification reactions was detected with the amount of template DNA used. Standard curves for each assay were generated by serial dilutions of linearized plasmids with cloned fragments of the respective gene. Standard curves were linear (R 2 = 0.997 ± 0.003) in the range used, and amplification efficiency was 90% for the 16S rRNA gene and 65–88% for the functional genes (**Supplementary Table S3**). Melting curve profiles were inspected, and final products were run on an agarose gel to confirm amplicon size. Non-template controls resulted in negligible values.

### Sequencing of the 16S rRNA Gene, Sequence Processing and OTU Clustering

The diversity and structure of total bacterial and archaeal communities were determined by targeting the V3-V4 region of the 16S rRNA gene (Takahashi et al., 2014). Amplicon libraries for each sample were generated using a two-step protocol (Berry et al., 2011). First, PCR products were generated in duplicate 20 µL reactions per sample using 16S rRNA primer constructs that included Nextera adapter sequences, with reactions consisting of Phusion PCR mastermix (Thermo Fisher Scientific), 0.5 µg µL −1 BSA, 0.25 µM of each primer and 10 ng extracted DNA. Thermal cycling was performed for 25 cycles, and cycling conditions and primer sequences are listed in **Supplementary Table S3**. The resulting PCR products were pooled and purified using the AMPure bead purification kit, and 3 µL of the purified product was used as template in the second PCR using barcodes. Duplicate 30 µL reactions were performed for each sample, with similar reagent concentrations as in the first step except for the use of 0.2 µM final primer concentrations. PCR was performed according to **Supplementary Table S3**. Products were pooled, bead purified, followed by equimolar pooling and sequencing performed by Microsynth AG (Balgach, Sweden) using the Illumina MiSeq platform with v3 chemistry (2 × 300 bp paired-end reads).

Paired-end reads were merged using PEAR (Zhang et al., 2013) and dereplicated and clustered into operational taxonomic units (OTUs) at a cut-off of 3% identity using UPARSE (Edgar, 2013). The final dataset comprised 13069 OTUs after removal of chimaeras and singletons, with 83% of the quality filtered sequence pool mapped back to OTUs. Taxonomic assignment

fmicb-10-01229 June 1, 2019 Time: 10:30 # 4

was carried out with the RDP classifier (Wang et al., 2007) against the SILVA reference database (release 119) (Quast et al., 2013). Sequences classified as mitochondria or chloroplasts were excluded. The original OTU table was rarefied (100 random subsampling) to 10660 sequences per sample. The sequences are available in the NCBI Sequence Read Archive (PRJNA494630).

### Statistical Methods

fmicb-10-01229 June 1, 2019 Time: 10:30 # 5

All multivariate, clustering and correlation analyses were performed using R (R Core Team, 2017). Comparisons of gene abundances between habitat types were performed using Kruskal–Wallis (KW) and Wilcoxon–Mann–Whitney (WMW) tests. Principal component analysis (PCA) and Redundancy analysis (RDA) using Hellinger distances (Borcard et al., 2011) were used to investigate the unconstrained ordination of the relative abundances of the N-functional genes studied (PCA) and of the bacterial community composition (PCA), and the relationship between the relative abundance of the N-functional genes and the environmental conditions (RDA), as well as between the overall bacterial community composition and the relative abundance of the N-functional genes (RDA). Hereafter, we refer to the PCAs as gene-PCA (gen-PCA) and community-PCA (com-PCA), and to the RDAs as gene-environment-RDA (gen-env-RDA) and community-gene-RDA (com-gen-RDA), respectively. In the analyses, functional gene abundances were standardised to total 16S rRNA gene copy numbers. Taxa < 5% occurrence (3453 of the total 13069 OTUs) were excluded from the bacterial community composition analysis (com-PCA and com-gen-RDA), and values for nitrate, nitrite, ammonium and sulphate were log-transformed in the gen-env-RDA. In the RDAs, forward selection was used to identify a minimum set of significant explanatory variables (p < 0.05; Blanchet et al., 2008), which exhibited low collinearity (variance inflation factors well below 10). Permutation tests of the resulting ordinations showed significant pseudo-F values (p < 0.05, n = 1000) for the main explanatory axes in each ordination — first to third axes in the gen-env-RDA, and first to fifth axes in the com-gen-RDA.

A structure of four sample clusters was present in both RDAs. Consequently, we used the samples scores of the three main axes of the com-gen-RDA as coordinates in four-group k-means clustering. We looked for indicative OTUs of each cluster performing a multi-level pattern analysis using the multipatt function from the indicspecies R package (Cáceres and Legendre, 2009), considering site group combinations, and the entire OTU set (13069) as the community data table. For each OTU, the method provides an indicator value (IndVal) of each cluster or a joint set of them. We accepted as significant indicator taxa those with adjusted p-value < 0.001, using the false discovery rate method to calculate the adjusted p-value (Storey et al., 2004).

### RESULTS

### Genetic Potentials

The sum of the N-functional gene copies per dry weight of sediment increased with OM (r = 0.42, p < 0.001, **Figure 2A**) but with a large scattering. Individual gene abundances were highly positively correlated among them (**Figure 2B**). However, the correlation structure markedly simplified when standardising by the 16S rRNA copy number in each sample (**Figure 2C**), showing that the nrfA pool was weekly related to the rest of N-functional gene pools. The sum of the N-functional gene copies accounted for an average of 15 ± 8% (mean ± SD) of bacterial 16S rRNA gene copies across all samples (**Supplementary Figure S1**). Maximum values of 52% were found in the lower sediment layer (2–4 cm) near the isoetid rhizosphere, while minimum values of 2% were observed in lithic biofilms. Hereafter, unless otherwise indicated, we report the N-functional gene abundance standardised to total bacterial 16S rRNA gene copies.

The nirS and nrfA genes showed the highest relative abundances, up to 33 and 18% of total 16S gene copies, respectively. The abundance of nirS was approximately 50-fold greater than nirK across all lakes, with the highest numbers detected in the more productive lakes (R, P and G, **Table 2** and **Figures 3A,B**). The abundance of nirK genes exhibited an overall trend of increasing abundance with lake DIN levels (**Figure 3B**), opposed to that observed for nirS. Higher nrfA abundance was observed in the sediments of the deepest part of the lakes (**Figure 3E**), while abundances in the elodeid sediments were significantly higher than those of the isoetids (KW test, p = 0.028). Closer inspection of the macrophyte sediment profiles showed a significant nrfA increase deeper in the elodeid sediments (r = 0.85, p < 0.001; **Supplementary Figure S2A**), while no trend was observed in the isoetid sediments.

The amoA gene of ammonia-oxidising archaea (AOA) was more abundant than the bacterial (AOB) counterpart. No obvious trend was observed for either AOA or AOB abundance across the lake DIN gradient. Although the average total abundance of ammonia oxidisers across all lakes was low relative to those of 16S rRNA genes (0.87 ± 2.66% and 0.05 ± 0.14% for AOA and AOB, respectively), several lakes showed AOA and AOB proportions of 57 and 23% of the total N-functional gene abundance. The highest AOA abundance was observed in the lower sediment layers of the isoetid rhizosphere and rocky littoral sediments of high-altitude lakes in the alpine belt, that is, those located above treeline (**Figure 3G**). Abundances of AOA were significantly higher in isoetid than elodeid sediments (KW test, p < 0.001), and increased with depth in the former (r = 0.64, p = 0.001; **Supplementary Figure S2B**). All habitats of the highest altitude lakes (GB and C lakes, **Table 2** and **Figure 3H**) showed a relatively high abundance of AOB copies compared to the same habitats in lakes at lower elevations.

The gene variants of the nitrous oxide reductase, nosZ clade I and II, as well as the hdh gene associated with the anammox pathway, exhibited low relative abundances (**Figures 3C,D,F**). The abundance of clade I nosZ genes was typically higher (∼7 fold on average) than that of clade II nosZ across all lakes and habitats. The lithic biofilm habitats of Contraix, the most elevated lake, showed the highest abundance of nosZ clade I genes (**Figure 3C**), while nosZ clade II abundances were higher in the rocky littoral sediments of alpine lakes (**Figure 3D**). Relative abundances of hdh were higher in the lithic biofilms (**Figure 3F**), with no obvious relationship with DIN levels across the lakes.

# N-Functional Genes and the Environment

The constrained ordination of N-functional gene abundances identified three distinct gradients explaining 54% of the total variation across habitats and lakes (**Figure 4**). A similar result was obtained in a non-constrained analysis (**Supplementary Figure S3**), indicating that the main environmental drivers were captured by the constrained analysis. The main variation of benthic N-cycling genetic potentials was across a nirS vs nrfA abundance gradient (**Figure 4**). The nirS-rich samples corresponded to those from shallow and productive lakes (R, P, and G, **Table 2** and **Supplementary Figures S4A,B**), specifically the upper sediments in all habitats and sediments near the roots of isoetids. These sites were associated with higher temperature, DNA content, isotopic signatures, DOC content, and C and N content, as well as coarser granulometry (**Supplementary Figure S5**) and lower C/N and nitrate/nitrite ratios (**Figure 4**). In contrast, nrfA rich environments occurred in the deep parts of the deep lakes, the lower layers of all sediments (except the isoetid rhizosphere) and the littoral sediments of the montane belt lakes (**Table 2** and **Supplementary Figures S4A,B**). Sulphate, ammonium, nitrate and nitrite concentrations were higher in these sites compared to those associated with nirS. The same nirS-nrfA main axis was also found if only samples from the deep habitat were included in the analysis.

The gen-env-RDA second axis of variation discriminated between sediments and lithic biofilms. The latter characterised by higher abundances of nirK, nosZI, nosZII, hdh and AOB (**Figure 4A**). The rocky littoral sediments of Lake Contraix separated from the other sediment samples (**Supplementary Figure S4A**). These sites shared relatively high concentrations of nitrate, nitrite and sulphate in the overlying water, high OM content, and particular isotopic signatures (high δ 13C and low δ 15N). Finally, the third axis of variation was associated with the AOA abundance (**Figure 4B**) and resulted in the segregation of the majority of the isoetid sediments from the other habitats. The sites with the highest abundance of AOA were located close to the isoetid rhizosphere and in the rocky littoral sediments of the alpine lakes (**Supplementary Figure S4B** and **Table 2**). These samples showed high δ 15N values and likely corresponded to more oxygenated sediments (**Figure 4B**).

# N-Functional Genes and the Associated Microbial Community

The ordination of the OTU composition constrained by the N-functional gene abundances resulted in a pattern of four distinct sample clusters (**Figure 5**), similar to that obtained in an unconstrained ordination (**Supplementary Figure S6**). The four clusters consisted of samples associated with a high relative abundance of nrfA, nirS, AOA, or a combination of the rest of

FIGURE 3 | Abundances of nirS (A), nirK (B), nosZI (C), nosZII (D), nrfA (E), hdh (F), archaeal amoA (G), and bacterial amoA (H) standardised to the 16S rRNA gene copies in different lakes and habitats. Numbers in brackets in the X-axis indicated the samples analysed. Boxplots depict the interquartile range (box), median value (line), 1.5x interquartile (whiskers), and outliers (points). Lakes are arranged from left to right by increasing water-column DIN level (see Table 2 for lake abbreviations and characteristics).

the targeted genes. Classification of samples into the four clusters using k-means followed by indicator species analysis resulted in approximately 29% of OTUs being identified as exclusively associated with samples from a single cluster (**Figure 6** and **Supplementary Tables S4, S5**). Approximately 12 % of OTUs were significant indicators of the AOA sample cluster and found

across a wide range of different bacterial taxa. By contrast, 6% of OTUs were significant indicators of the mixed N-transformation cluster, with large numbers of indicators concentrated within the phyla Cyanobacteria, Bacteroidetes, and Planctomycetes, as well as Alpha- and Betaprotoebacteria classes. Similarly, 6% of OTUs were associated with samples in the nrfA cluster and were classified into Firmicutes, Bacteroidetes, Actinobacteria, and Chloroflexi phyla, and Epsilon- and Deltaproteobacterial classes. Finally, 4% of OTUs were exclusive indicators of samples in the nirS cluster and were found across a large number of bacterial taxa, similar to the AOA sample cluster.

# DISCUSSION

### The DNRA-Denitrification Gradient

The gradient of nrfA to nirS dominance was the main pattern of variation in the N-transforming microbial communities of the benthic habitats. From an ecosystem perspective, this gradient indicates a shift from habitats with a higher potential for internal Nr cycling via DNRA, and thus retention of Nr in the system, to those in which loss of Nr from the lake is more likely through denitrification. The environments with the lowest ratios of denitrifying to DNRA nitrite reductase genes were characterised by variables indicating refractory OM with high C/N, lower oxygen diffusion and lower redox potentials. This was particularly the case in the deepest part of the deep lakes (maximum depth ≥ 25 m) and in the deeper regions of the reduced elodeid sediments. In agreement with our result, previous work in a tropical high-altitude oligotrophic lake has shown nrfA abundance to be highest in the deepest part of the hypolimnion with anoxic conditions during late stratification (Pajares et al., 2017), and more reduced conditions favoured DNRA over denitrification in Australian estuaries (Kessler et al., 2018). Increased C/N ratio may favour DNRA over denitrification (Kraft et al., 2014).

The highest genetic potential for denitrification, based on nirS abundance, was detected in shallower and less oligotrophic lakes, where DIN levels were lower in the water column due to higher primary productivity. These lakes also have lower C/N ratios (**Supplementary Table S1**). Generally, lake autochthonous OM is fresher and of higher quality (e.g., lower C/N), and is a substantial proportion of total OM in lakes with a lower ratio of the catchment to the lake area. This fresh OM can be used as electron donors for denitrifiers, as demonstrated in several aquatic ecosystems [eutrophic lakes (Chen et al., 2012; Gardner et al., 2017), streams (Barnes et al., 2012; Stelzer et al., 2014), wetlands (Dodla et al., 2008), and oceans (Van Mooy et al., 2002)]. Oxygen levels in upper sediments of shallow and productive lakes likely fluctuate to a greater degree than those observed in habitats dominated by DNRA, thereby favouring organisms with facultative anaerobic respiration pathways such as denitrification (Wittorf et al., 2016; Chen J. et al., 2017). The nirS denitrifiers along the DNRA-denitrification gradient were associated with nosZ clade II N2O reduction and AOA communities involved in ammonia oxidation, whereas denitrifier communities present in lithic biofilms dominated by nirK-types were associated with nosZ clade I N2O reduction and AOB. These patterns indicate

that different N transformation networks developed in these habitats even when in both cases exhibited potential for linked nitrification and denitrification.

The overall high proportion of denitrifying nitrite versus nitrous oxide reductase genes (∼30 nir:nosZ ratio on average) suggests a dominance of partial denitrification, especially in productive habitats dominated by nirS denitrifier communities (i.e., nirS-cluster showed nir:nosZ higher ratios compared to the other clusters, WMW test, p < 0.01). This observation agrees with Castellano-Hinojosa et al. (2017) who found high N2O/N2 emissions in a productive, shallow warm Mediterranean mountain lake, as well as Myrstener et al. (2016) who demonstrated that addition of nitrate, phosphorus and labile C to sediments from a boreal lake resulted in higher relative N2O production compared to addition of nitrate alone. Other studies have shown that higher nir:nosZ1 ratios in the sediments of boreal lakes were associated with hypolimnion N2O excess, as well as increased phosphate and nitrate concentrations (Saarenheimo et al., 2015a). Thus, productive sites could favour partial denitrifiers that survive anoxic periods. Before arriving to the atmosphere, N2O might be consumed in the hypolimnion of deep lakes. NosZ-harbouring bacteria have been found in the hypolimnion of boreal lakes (Peura et al., 2018). However, in the sediments studied, the higher nir:nosZ ratios were found in shallower lakes. in which N2O may easily reach the atmosphere (Dore et al., 1998). Further studies accounting for real N2O emissions could corroborate our conjecture.

Nitrite-dependent anaerobic methane oxidation (N-DAMO, Simon and Klotz, 2013) is a potential alternative to DNRA, denitrification and anammox nitrite consumption. N-DAMO

has been found as a key driver of methane oxidation in nitrate-rich lakes (Deutzmann et al., 2014) and reduced sandy riverbeds (Shen et al., 2019). In these habitats, bacteria related to Candidatus Methylomirabilis oxyfera, known to perform this pathway, were abundant. In our study, we detected only two OTUs of low relative abundance (0.02%) classified as Candidatus Methylomirabilis. Nonetheless,

clusters, as delimited by the dominant N-functional gene abundances.

OTUs belonging to the genus Methylocaldum were more abundant (∼0.3%), and were significant indicators of the 'nirS' cluster. These findings are similar to those of a survey of methane oxidation in Indian reservoirs (Naqvi et al., 2018), where low relative abundance of NC10 bacteria capable of N-DAMO (0.003–0.022%) was found; whereas Type I aerobic methanotrophs, which include Methylocaldum and other

members of the Methylococcaceae family, were predominant and co-occurred with a diverse community of potential nirS type denitrifiers. Other Methylococcaceae are partial denitrifying aerobic methanotrophs with N2O as the end product (Kits et al., 2015a,b). Overall, these results, added to the high ratio of nir to nosZ gene abundance, suggest that N2O emissions are the most likely endpoint of nitrite reduction in nirS-cluster sediments, independent of the pathway. Further studies are required to elucidate the importance and distribution of methane dependent processes in mountain lakes and to evaluate their role as a bypass of partial denitrification.

### The Idiosyncratic Lithic Biofilms

Microbial N-transforming guilds in the lithic biofilms differentiate from those in the sediments. Gene abundance results indicate a complex N-transformation structure in this habitat, consisting of processes demanding both oxic and anoxic conditions, which suggest highly structured microbial communities in relatively short spatial distances. The idiosyncratic nature of the guild composition (e.g., nosZI, hdh, nirK) declines as the productivity of the lake increases; nirS becomes more prominent compared to nirK, and the N-transforming communities are more similar to the upper sediments. The high abundance of nirS differs from the dominance of nirK previously found in another study of epilithic biofilms from a subset of the same lakes. A main difference between the two studies is that in Vila-Costa et al. (2014) only sampled the upper side of the cobbles, whereas we sampled both sides. Differences in ammonia oxidisers between the upper (light-side) and lower (dark-side) sides of cobbles have been previously reported (Merbt et al., 2017). The higher relative abundance of hdh and nosZ genes indicates that N loss in the form of N2 could be higher in the lithic biofilms compared to the other benthic habitats in the studied lakes.

### Archaeal Nitrification Hotspots

The positive δ 15N signals observed in samples from the lower sediment layers of the isoetid rhizosphere and the rocky littoral sediments of the alpine lakes support the view of them as nitrification hotspots (Mariotti et al., 1981), likely performed by AOA as suggested by the amoA gene abundance. Nitrification in rocky littoral sediments could be quantitatively more relevant than the isoetid rhizosphere given that rocky littorals occupy large areas in alpine lakes. Nevertheless, the archaeal amoA densities observed in the lower isoetid sediment layer were nearly 100-fold higher than in the rocky sediment samples with the highest AOA abundance. Sediments close to the roots of isoetids are episodically well-oxygenated due to the release of oxygen through roots during photosynthetic periods (Sand-Jensen et al., 1982), which increases the interface between oxidised and reduced sediments where the NH4 + oxidation occurs. Intensive nitrification would result in high NO3 − accumulation in the sediment porewater (maxinun 2 mM; Catalan et al., 1994). However, this is likely transient as the NO3 − concentration in the overlying water column was negligible, suggesting a close coupling of nitrification and denitrification in this habitat (Vila-Costa et al., 2016). Indeed, there was a positive correlation (r = 0.68, p < 0.001) between nitrification and denitrification gene abundances in the isoetid sediments of Plan lake. Lakes at a higher altitude, Contraix and Gelats de Bergús, also showed significant correlations between nitrification and denitrification gene abundances for deep and littoral habitats (r = 0.48, p = 0.02 and r = 0.44, p = 0.08, respectively), suggesting that nitrification and denitrification are also linked in the nitrification hotspots.

# Linking the Taxonomic Distribution and Functional Potential

Each of the four main N-transforming communities consists of a highly diverse and distinct consortium of co-occurring bacterial taxa, based on the large number of indicator OTUs detected. However, taxonomic classification does not necessarily predict functioning, as different prokaryotic traits may be conserved at different phylogenetic depths (Martiny et al., 2015). While ammonia oxidation and anammox capacities are restricted to only a few lineages, denitrification and DNRA are widely distributed across the phylogeny (Graf et al., 2014; Welsh et al., 2014). Thus, many of the indicator OTUs identified are likely not directly involved in each N-transformation pathway. However, the high degree of similarity observed between unconstrained and functional gene-constrained OTUbased ordinations indicates that shifts in the genetic potential of different N-transformation processes are tightly linked to changes in the overall prokaryotic community structure, which itself is shaped by differences in environmental conditions across habitats.

Links between taxonomic composition and functional potential has been observed in previous works in lakes based on metagenomes or sequencing of functional genes. The occurrence of several proteobacterial families, in particular, Rhodobacteraceae (Alphaproteobacteria), Methylococcaceae (Gammaproteobacteria), and Burkholderiales, Comamonadaceae and Rhodocyclaceae (Betaproteobacteria) has been shown to be strongly associated with denitrification gene presence or abundance (Vila-Costa et al., 2014; Peura et al., 2015; Saarenheimo et al., 2015b; Castellano-Hinojosa et al., 2017; Chen R. et al., 2017). These taxa were also highly abundant in samples within the nirS-denitrifier and mixed functional gene communities, which included nirK-type denitrifiers. Moreover, metagenomic studies of boreal lakes water columns have identified nosZ sequences originating from Myxococcales (Deltaproteobacteria) and Sphingobacteriaceae (Bacteroidetes) in the hypolimnion near the oxycline (Peura et al., 2015, 2018). Many organisms within these families are known to possess the clade II variant of the nosZ gene (Hallin et al., 2018). Accordingly, a large proportion of indicator OTUs for the mixed N-cycling communities, which includes the clade II nosZ variant, were also classified as belonging to these families. The 4th most abundant genus in the studied lakes is Anaeromyxobacter, one member of this Deltaproteobacteria genus is A. dehalogenans a chemodenitrifier, an organism that combines chemical chemodenitrification reactions and enzymatic reaction(s) to reduce NO3 − to N2O or N2, without

having denitrifying nitrite reductases codified by nirS or nirK, also performs DNRA and Fe-reduction (Onley et al., 2018). Rhodoferax (Beta-) and Desulfomonile (Delta-) possible chemodenitrifiers (Onley et al., 2018) were also common genus present. There are other eubacteria non-proteobacteria taxa also carrying denitrifying genes (Graf et al., 2014) present in our samples (e.g., Actinobacteria).

Operational taxonomic units associated with samples in the nrfA cluster were classified as Firmicutes, Epsilon- and Deltaproteobacteria (Campylobacterales and Anaeromyxobacter, Desulfovibrio, and Geobacter, respectively), Bacteroidetes (Bacteroidia), Actinobacteria (Coriobacterales and Corynebacterales) and Chloroflexi (Anaerolineaceae), all these taxa include microbes that are known to carry nrfA (Welsh et al., 2014).

Regarding the ammonia oxidisers, the primers used in the 16S rRNA sequencing mainly target bacteria, but also pick up Euryarchaeota (Takahashi et al., 2014). Therefore, OTUs assigned to Thaumarchaeota were not detected, although AOA hotspots could be identified based on qPCR data. Nitrosomonadaceae was the most common AOB, with 89 OTUs classified as being similar to uncultured members of this family. Nitrospirae was the only identified nitrite oxidising bacteria (NOB) in our samples. There was a likely coupling between AOA and NOB, as suggested by the correlation between archaeal amoA genes and the relative abundance of Nitrospirae members in general, and Nitrospira in particular, in the AOA-cluster samples (r = 0.65, p = 0.0005; r = 0.47, p = 0.02, respectively). This coupling has previously also been found in grasslands (Simonin et al., 2015), agricultural soils (Jones and Hallin, 2019), and sediments of an Andean mountain lake (Parro et al., 2019). Comammox Nitrospira could be important in the nitrifying hotspots found in the present study, as suggested by previous studies in other surface-attached oligotrophic habitats (Kits et al., 2017; Pjevac et al., 2017; Fowler et al., 2018). For anaerobic ammonia oxidation, OTUs belonging to the "Candidatus Anammoximicrobium" (Khramenkov et al., 2013) was the only taxa present with demonstrated anammox capacity. However, more bacteria within the numerous uncultured Planctomycetes detected in our samples could potentially perform anammox.

### CONCLUSION

The N-transforming guild composition in benthic habitats of mountain lakes is complex and deeply embedded in the overall prokaryotic community. There is a high positive correlation among all the genes, and they all generally increase with OM. The dominant pathways change depending on the habitat and productivity of the lake (**Figure 7**). The fate of nitrite is the main diverging point differentiating the N-transforming guilds. The genetic potential for DNRA dominate in the deep part of the lakes and the lower sediment layers, which indicates recycling of the Nr . By contrast, the denitrifying nirS nitrite reduction potential prevails in the upper layer of the sediments in the shallow, warmer and more productive lakes, which indicates a loss of Nr . Emissions of N2O and N2 are likely spatially segregated within lakes, with lithic biofilms being candidates for preferential N loss as N2 as they show a more balanced gene abundance of nitrous oxide reductases (nosZI+II) and anammox (hdh) in relation to NO-forming nitrite reductases (nirS+nirK). The more productive and nirS-dominated habitats may be a main source of N2O because of the striking excess of this gene over the ones of the final steps of complete denitrification unless another bypass process is relevant (e.g., N-DAMO). There may be two types of nitrifying-denitrifying coupled community types in the benthic habitats of mountain lakes. The first is based on nitrification by AOA coupled to Nitrospirae (NOB) and denitrification by nirS-denitrifiers, with hotspots in the rocky littoral sediments of the lakes above treeline and the sediments near the isoetid rhizosphere. The second includes AOB coupled to nirK-type denitrifiers reducing nitrite and nosZI-N2O reduction in the lithic biofilms. Overall, our results point out two types of potential response to high atmospheric N deposition in these lakes. In highly oligotrophic lakes, there will be an accumulation of Nr because of the predominance of internal Nr recycling via DNRA. In less oligotrophic lakes, generally with macrophyte growth, the Nr deposition loads may be more effectively directed toward N gas release to the atmosphere via denitrification.

### ETHICS STATEMENT

fmicb-10-01229 June 1, 2019 Time: 10:30 # 13

The authors declare that the present study does not involve human or animals, that they have all the licenses, and that they follow all the rules for sampling in the Aigüestortes National Park.

### AUTHOR CONTRIBUTIONS

CP-L and JCat contributed to the study design. LC and CP-L carried out sampling. CP-L, CJ, and JCal carried out the lab work and data analysis. JCat, SH, LC, and EC contributed to reagents, materials, and analysis tools. CP-L and JCat wrote the manuscript. All authors substantially contributed to commenting and revising it.

### REFERENCES

### FUNDING

The Spanish Government provided funds through the Ministerio de Educación as a predoctoral fellowship to CP-L (FPU12–00644) and research grants of the Ministerio de Economia y Competitividad: NitroPir (CGL2010–19737), DARKNESS (CGL2012-32747), LACUS (CGL2013–45348-P), TRANSFER (CGL2016–80124-C2-1-P). The REPLIM grant (INRE - INTERREG Programme. EUUN – European Union. EFA056/15) and the ALTER-NET Mobility Scheme supported the final writing.

### ACKNOWLEDGMENTS

We thank the authorities of the Aigüestortes and Estany de Sant Maurici National Park for sampling facilities in protected areas and support.

### SUPPLEMENTARY MATERIAL

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

of bacterial nitrate respiration. Science 345, 676–679. doi: 10.1126/science. 1254070

denitrifying bacteria with DGGE. FEMS Microbiol. Ecol. 49, 401–417. doi: 10. 1016/j.femsec.2004.04.011

**Conflict of Interest Statement:** 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.

Copyright © 2019 Palacin-Lizarbe, Camarero, Hallin, Jones, Cáliz, Casamayor and Catalan. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Decoupling the Dynamics of Bacterial Taxonomy and Antibiotic Resistance Function in a Subtropical Urban Reservoir as Revealed by High-Frequency Sampling

Peiju Fang1,2, Feng Peng1,2, Xiaofei Gao1,2, Peng Xiao1 and Jun Yang1 \*

1 Aquatic Ecohealth Group, Key Laboratory of Urban Environment and Health, Institute of Urban Environment, Chinese Academy of Sciences, Xiamen, China, 2 University of Chinese Academy of Sciences, Beijing, China

### Edited by:

Haihan Zhang, Xi'an University of Architecture and Technology, China

### Reviewed by:

Wenguang Xiong, South China Agricultural University, China Yihui Ban, Wuhan University of Technology, China

> \*Correspondence: Jun Yang<EMAIL_ADDRESS>
### Specialty section:

This article was submitted to Aquatic Microbiology, a section of the journal Frontiers in Microbiology

Received: 03 April 2019 Accepted: 11 June 2019 Published: 02 July 2019

### Citation:

Fang P, Peng F, Gao X, Xiao P and Yang J (2019) Decoupling the Dynamics of Bacterial Taxonomy and Antibiotic Resistance Function in a Subtropical Urban Reservoir as Revealed by High-Frequency Sampling. Front. Microbiol. 10:1448. doi: 10.3389/fmicb.2019.01448 Aquatic environments serve as important reservoirs of antibiotic resistance genes (ARGs), but the information on the high-resolution temporal pattern of ARGs in waterbodies is extremely limited. In this study, the weekly dynamics of ARGs and their relationships with microbial taxonomic communities and environmental variables were analyzed in a subtropical urban reservoir over the period of 1 year using highthroughput approaches. In total, 197 ARGs and 10 mobile genetic elements (MGEs) were detected. The results showed that the bacterial community had a seasonal pattern, while ARGs composition did not exhibit seasonality, thereby indicating the asynchrony or decoupling of temporal patterns of microbial taxonomy and function. More importantly, bacterial abundance and community diversity were more strongly correlated with 17 measured environmental variables than ARGs (36 significant correlations for OTUs, 11 for ARGs). However, stochastic processes appeared to have a minor role in the structuring of the ARG profiles, but a more important role in the structuring of bacterial taxonomic communities. Furthermore, we found that precipitation and turbidity were significantly correlated with the richness and diversity of ARGs, suggesting that multiple environmental factors influence the composition and dynamics of ARGs in complex ways. MGEs were abundant and showed significant positive correlations with ARGs, indicating a plausible influence of MGEs on the variation of ARGs. This is the first study which provides an overview of high-resolution dynamics of ARGs in a subtropical waterbody. Our results improve the understanding of microbial processes and mechanisms of ARGs at fine temporal scale, and offer empirical data of use in the monitoring, assessment and management of the urban water environments.

Keywords: antibiotic resistance genes, bacterioplankton, microbial community, high-resolution temporal pattern, subtropical reservoir

# INTRODUCTION

In aquatic ecosystems many studies have been carried out on nitrogen, phosphorus, and microbial community (Yang et al., 2012; Zhang et al., 2019), however various "emerging" pollutants including antibiotics and antibiotic resistance genes (ARGs) have been rapidly drawing increased attention from the public and government (Pruden et al., 2006; Zhu et al., 2017; Yang et al., 2018).

Antibiotics, first discovered by Sir Alexander Fleming in 1928, were regarded as panacea for any microbial infections (World Health Organization [WHO], 2014). However, with the emerging problem of antibiotic resistance, antibiotic resistant bacteria (ARB) and ARGs have been identified as a significant issue, posing a serious threat to the health of humans and the environment (Pruden et al., 2006; Berendonk et al., 2015). Antibiotic resistance is ancient (D'Costa et al., 2011), but with the continuous increase in discharge of overused antibiotics, the development, spread, and enrichment of ARGs has accelerated globally (Laxminarayan, 2014), and heavy metal has been reported to play an important role in the enrichment of antibiotic resistance by co-selection mechanisms (Stepanauskas et al., 2006). More importantly, ARGs can be exchanged and transferred among environmental bacteria, pathogens and nonpathogens via horizontal gene transfer (HGT) carried out by mobile genetic elements (MGEs) (Kruse and Sørum, 1994; Martínez et al., 2015). Aquatic environments are considered to be significant reservoirs of both ARB and ARGs, providing ideal settings for the occurrence and dissemination of ARGs (Marti et al., 2014). Due to long-term misuse and overuse of antibiotics by human, urban water often contains diverse ARGs at high levels (Xu et al., 2016), which may pose high risk to human health during the interaction with receivers. Therefore, ARGs occurrence and their temporal dynamics in the aquatic environment have become an important issue in environmental science (Yang et al., 2018).

Microbial communities can vary between environments, and play an important role in drive global biogeochemical cycling, but it is not well known how variation in taxonomic composition relates to function (Louca et al., 2016). Previous studies reported that bacterial community composition closely correlates with ARG profiles and controls the transfer of ARGs in antibiotic-rich environments (e.g., soil, urban wastewater, and sewage sludge) (Forsberg et al., 2014; Su et al., 2015), but the contribution of microbial community shifts to antibiotic resistance function variations in natural waterbodies with low and medium antibiotic concentrations remains largely unknown, especially at a fine temporal scale. Furthermore, microbial communities have been reported to exhibit seasonal variability in the environment (Treusch et al., 2009), and high-frequency sampling analysis could reveal and disentangle detailed temporal dynamics of bacterial community (Lindh et al., 2015). Previous studies have assessed the temporal variation of ARGs and supported the occurrence of seasonal variations of ARGs based on only two or four sampling time points in 1 year (An et al., 2018; Zheng et al., 2018). However, research on the dynamics of ARGs in urban water environment, and their relationship to microbial taxonomic community and environmental factors based on highfrequency time series is still lacking.

Revealing the mechanisms that drive the variation of microbial communities and ARG profiles are major challenges. It has been shown that microbial communities are simultaneously shaped by stochastic (neutral) and deterministic processes (Sloan et al., 2006; Logares et al., 2013). The neutral model is a stochastic process that develops community through births, death and immigration, can successfully predict microbial community assembly dynamics (Hubbell, 2001). However, the neutral model cannot completely descript the community assembly without deterministic process such as environmental factors and competition (Sloan et al., 2006; Ofi¸teru et al., 2010). Both stochastic (neutral) and deterministic processes need to be considered when investigating the assembly of microbial communities and ARG profiles, but there is still a lack of studies on the contribution of both processes in structuring ARG profiles (Zhu et al., 2017; Guo et al., 2018).

The urban reservoir in this study is located in a rapidly urbanizing area of Xiamen, it receives numerous pollutants from surrounding housing estates by point and non-point source pollutions. However, this urban reservoir is the most important landscape waterbody in the Jimei district, so the ARB and ARGs in the water may pose ecological risk in the interaction with inhabitants and tourists. In this study, we simultaneously investigated the dynamics of bacterial OTUs and ARGs in the urban reservoir using high-throughput qPCR and 16S rRNA gene sequencing on a weekly basis over a 1-year period. We aimed to (1) characterize the temporal dynamics of both bacterial taxonomic community and antibiotic resistance gene profiles; (2) reveal the relationships between bacterial taxonomic composition and antibiotic resistance function at different levels of taxonomic resolution; and (3) explore the community assembly processes and mechanisms shaping the ARGs composition dynamics. Additionally, we hypothesize that the temporal patterns of microbial OTUs and ARGs are asynchronous, and their temporal patterns are influenced by different processes and mechanisms. To our best knowledge, this is the first study to analyze a high-resolution time series of bacterial OTUs and ARGs communities in urban water environments. Therefore, this study may provide a fundamental data for evaluating the ecological risks of ARGs in urban water environment and offer useful information for government and policy maker.

# MATERIALS AND METHODS

### Study Area and Sample Collection

All samples were collected once a week from September 2016 to August 2017 in station G (24◦ 360 N, 118◦ 040 E) of Xinglinwan Reservoir, which is situated in the Houxi River watershed, Xiamen city, Fujian province, southeast China (**Supplementary Figure S1**). In total, 51 samples were collected in this study. Several stormwater and wastewater inlets were found near the station G along the reservoir shoreline within a distance of 1000 m. The Xinglinwan Reservoir is a subtropical eutrophic urban reservoir, and is characterized by high turbidity and low water transparency. It has multiple functions in flood control, aquaculture and tourism, exerting a great impact on the surrounding residents.

Surface waters were collected in a sterile container from the top 0.5 m water, immediately transported to laboratory and kept in the dark until pre-treatment within 0.5 h. In order to remove the influence of large plankton and particles, samples were prefiltered through a 200 µm mesh, and then filtered through

0.22 µm polycarbonate membranes (47 mm diameter, Millipore, Billerica, MA, United States) using a vacuum filtration system. To ensure the sufficient microorganisms, the filtrated water volume ranged from 200 to 500 mL for all samples, because we kept a similar filtrating time (about 50 min) for each of membranes. The membranes were placed in sterilized tubes and stored at –80◦C until DNA extraction.

### Environmental Parameters

A total of 17 environmental variables were determined in this study (**Supplementary Figure S2**). The comprehensive trophic state index (TSIc) were calculated according to Yang et al. (2012). Water temperature (WT), pH, dissolved oxygen (DO), turbidity, electrical conductivity (EC), salinity, and oxidation reduction potential (ORP) were measured in situ with a Hydrolab DS5 multiparameter water quality analyzer (Hach Company, Loveland, CO, United States). Chlorophyll a concentrations were quantified by PHYTO-PAM Phytoplankton Analyzer (Heinz Walz GmbH, Eichenring, Germany). Total carbon (TC), total organic carbon (TOC), total nitrogen (TN), ammonium nitrogen (NH4-N), nitrate nitrogen (NO3-N), nitrite nitrogen (NO2- N), total phosphorus (TP), and phosphate phosphorus (PO4-P) were measured according to the standard methods described in our previous study (Liu et al., 2013). The precipitation data referred to the cumulative rainfall during the last week before the sampling day, and the data were collected from the Xiamen Meteorological Bureau.

# DNA Extraction and High-Throughput Sequencing

Total DNA was directly extracted from each membrane using the FastDNA SPIN Kit (MP Biomedical, Santa Ana, CA, United States) according to the manufacturer's protocol, and subsequently the concentration and quality of the DNA were determined using a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States).

To investigate the bacterial taxonomic community, the V3- V4 hypervariable regions of the 16S rRNA gene were amplified, purified, and quantified according to our previous procedure (Guo et al., 2018). The PCR products were pooled and sequenced on the Illumina Miseq platform (Illumina, Inc., San Diego, CA, United States) using a paired-end (2 × 250 bp) sequencing strategy. All the raw sequences were quality-controlled using MOTHUR v1.39.0 (Schloss et al., 2009). The unoise3 pipeline was used to pick operational taxonomy units (OTUs) at 3% dissimilarity level (Edgar, 2010). OTU sequences were taxonomically classified by USEARCH (sintax) against the Greengenes database (DeSantis et al., 2006). All eukaryotic, chloroplast, archaeal, mitochondrial and unknown sequences were removed from the data set. Finally, all sequences were normalized to 50,014 sequences for each of 51 samples, and these sequences were clustered into 9834 OTUs.

The raw data for 16S rRNA gene analysis in this study have been deposited in the NCBI sequence read archive (SRA) database under the BioProject number PRJNA510463 and the accession number SRP173857.

# High-Throughput Quantitative PCR (HT-qPCR)

High-throughput quantitative PCR was performed to detect the abundance of ARGs using the Wafergen SmartChip realtime qPCR platform (Wafergen Biosystems, Fremont, CA, United States). A total of 296 primer sets targeting 285 resistance genes, 8 transposase genes, the class I integron-integrase gene (intI), the clinical class 1 integron-integrase gene (cIntI), and 16S rRNA gene were used. These primer sets and PCR reaction protocols were listed in our previous studies (Guo et al., 2018; Liu et al., 2018). Three technical replicates were performed for each sample, and non-template negative control was also included for each primer set. The absolute abundance of ARGs and MGEs was calculated by normalizing their abundance to the absolute 16S rRNA gene copy number, which was estimated by real-time quantitative PCR (Ouyang et al., 2015) as described below.

# Real-Time Quantitative PCR (qPCR)

Real-time quantitative PCR was used to quantify the absolute copy number of 16S rRNA genes in all samples. It was performed on a Lightcycler 480 instrument (Roche, Basel, Switzerland). The 20 µL reaction contained 10 µL 2× LightCycler 480 SYBR Green I Master Mix (Roche, Basel, Switzerland), 7 µL nucleasefree PCR-grade water, 2 µL diluted DNA, and 0.5 µM of each primer. The reactions were performed in triplicate with negative controls. The following thermal cycling conditions were used: initial incubation at 95◦C for 5 min, followed by 40 cycles of 95◦C for 15 s, 60◦C for 1 min, and 72◦C for 15 s.

A six-point calibration curve was generated from 10-fold dilutions for an external standard calculation. The melting curves were used to analyze the specificity of PCR products, and PCR efficiency ranged from 96 to 105% in this study.

# Statistics

The Bray-Curtis similarity was calculated using PRIMER v7.0 (Clarke and Gorley, 2015), and non-linear regression of time lag with community similarity was analyzed to reveal the temporal patterns of bacterial OTUs and ARGs communities. The nonlinear regression was the fit to the sinusoidal curve to explore the seasonality in SigmaPlot v12.0 (Systat Software Inc., Chicago, IL, United States). The diversity indexes were calculated using the vegan package in R v3.3.1. The SIMPER analysis was performed to identify the contribution of bacterial OTUs and ARGs to the temporal variation of communities or profiles (Clarke and Gorley, 2015).

Spearman's and Pearson's correlations were calculated using SPSS v20.0 (IBM Corp, Armonk, NY, United States), and the Mantel test was performed in an R environment with the vegan package. Further, a Procrustes test was conducted with the vegan package in R environment to explore the synchronicity of bacterial OTUs communities and ARG profiles (Forsberg et al., 2014). Network analysis was used to reveal Spearman correlations between ARGs and bacterial taxa with the picante package in R v3.3.1, and the correlation coefficients |ρ| > 0.6 between microbial taxa and ARGs, and |ρ| > 0.4 between MGEs and

ARGs, were considered at P < 0.01. Network information was generated by the Gephi v8.2 and Cytoscape v3.6.1.

In addition, we used redundancy analysis (RDA) to explore the relationship between environmental factors and the compositions of bacterial OTUs and ARGs subtypes. We further used partial redundancy analysis (pRDA) to evaluate the relative contribution of bacterial taxonomic community, environment variables and MGEs to the ARG profiles of all samples by R v3.3.1 with the vegan package. To reveal the assembly processes underlying both bacterial OTUs and ARGs compositions, we used the neutral community model, which can evaluate the importance of stochastic processes in explaining the assembly of bacterial taxonomic and functional communities (Sloan et al., 2006). This model predicts the relationship between the frequency at which taxa occur in a set of local communities (in this case, bacterial OTUs and ARGs at one sampling point) and their abundance across the wider metacommunity (bacterial OTUs and ARGs sampled over a year). The neutral models were performed in R v3.3.1 using the Hmisc, stats4, and minpack.lm packages.

## RESULTS

### Dynamics of Environmental Variables, Bacterial OTUs, and ARGs Compositions

The major physicochemical variables exhibited a great change over time (**Supplementary Figure S3**). Comprehensive TSIc varied from middle eutrophic level to hypereutrophic level, and nutrient concentrations were high throughout the study period, with higher value in 2017 than 2016. Water temperature ranged from 15 to 35◦C during the 1-year sampling. Higher precipitation values were recorded in September 2016 and June 2017.

We found 197 unique ARGs and 10 MGEs in this study, and the ARGs covered almost all major classes of antibiotics commonly used. The number of ARGs detected in each sample varied from 18 to 74, representing three major resistance mechanisms including antibiotic deactivation (53.8%), cellular protection (28.6%), and efflux pump (16.9%) (**Supplementary Figure S4**). The absolute abundance of ARGs ranged from 4.03 × 106 copies/L to 3.72 × 109 copies/L (mean ± s.e. = 4.55 × 108 ± 1.11 × 108 ), while the absolute abundance of MGEs varied from 6.62 × 105 copies/L to 1.84 × 109 copies/L (1.92 × 108 ± 5.28 × 107 ). The Shannon-Wiener index of bacteria ranged from 4.95 to 7.58 (**Supplementary Figure S5**). The bacterial community was dominated by Actinobacteria, Proteobacteria and Cyanobacteria, while Proteobacteria was the most abundant phylum in September and October (36.1% ± 2.5%, n = 8), and Actinobacteria dominated the bacterial community most of the time (37.9% ± 1.7%, n = 43).

Time series analysis indicated relative stability at bacterial phylum and ARG type level, but very rapid or abrupt fluctuations among both bacterial OTUs and ARG subtypes (**Figure 1**). The relative abundance of ARG subtypes exhibited distinct variations across the year at weekly intervals (**Figure 1D**). For example, the relative abundance of qacE11-02 (multidrug resistance gene) varied from 0 to 0.5 and its relative abundance changed almost every week. Time-lag regression analysis showed that the bacterial community displayed a good fit to a sinusoidal curve (**Figure 1E**), but ARG profiles failed to fit such a curve (**Figure 1F**); indicating that bacterial community composition exhibited seasonality while ARGs composition did not. The general variation of ARG similarity was narrower compared to the bacterial taxonomic community (**Supplementary Figure S6**), suggesting that the compositions of resistomes were more stable. Additionally, SIMPER analysis indicated that Proteobacteria contributed most to the variation within the microbial community (31.86%), while beta-lactamase resistance gene contributed most to the variation of ARGs (23.41%) (**Supplementary Table S1**). Furthermore, Procrustes test revealed asynchrony or decoupling between bacterial ARGs and OTUs communities (r = 0.347, P < 0.01, **Figure 2A**). However, the Pearson's correlation indicated a significant but weak relationship between bacterial ARGs and OTUs communities (**Figure 2B**).

### Co-occurrence Patterns Between Bacterial Taxonomic Community and Antibiotic Resistome

Network analysis was performed to reveal the significant and robust correlations between ARG subtypes and bacterial OTUs. The whole network consisted of 1,547 nodes and 2,139 edges, and showed a modularity index of 0.762. Aminoglycoside resistance genes had the most edges (489) with bacterial taxa, followed by tetracycline (474), MLSB (454), and beta-lactamase (318) resistance genes (**Figures 3A,B**). Meanwhile, Proteobacteria, Bacteroidetes, Firmicutes, and TM7 were significantly correlated with most ARGs. The co-occurrence network showed that the hubs of all modules were ARGs (**Supplementary Figure S7A**), indicating the important roles of ARGs in the network.

We also used Mantel tests to evaluate the relationships between bacterial phyla and ARG types. Similar to the results of the network analysis, most ARGs were significantly and strongly correlated with Bacteroidetes, Firmicutes, and Proteobacteria (**Supplementary Table S2**). However, sulfonamide resistance genes did not show any significant correlation with the bacterial community, and MGEs only significantly correlated with Chloroflexi, Firmicutes, and TM7.

### Relationships Between Environmental Factors, Bacterial OTUs, and ARGs

Mantel tests showed aminoglycoside, beta-lactamase, MLSB, tetracycline, and vancomycin resistance genes were significantly correlated with precipitation (**Supplementary Table S3**). RDA further indicated that water temperature, chlorophyll a, salinity, total carbon, TOC, total nitrogen, nitrate nitrogen, and total phosphorus were significantly correlated with the β-diversity of the bacterial community, while salinity, total carbon, nitrate nitrogen, total phosphorus, and precipitation were strongly correlated with ARGs compositions, suggesting that nutrients were key factors that influence the dynamics of

both ARGs and bacterial communities (**Figure 4**). We found that the bacterial community (R 2 = 0.365) showed a much better fitted to the neutral model (**Figure 4**) than ARGs (R 2 = 0.079), indicating that stochastic process played an important role in microbial community assembly but not in the resistome profiles.

The absolute abundance of ARGs was positively correlated with total carbon and TOC, while ARGs richness was significantly and positively correlated with turbidity and precipitation (**Figure 5A**). Additionally, turbidity and precipitation also exhibited a positive correlation with richness and α-diversity of bacterial community, but were not significantly correlated with 16S rRNA gene abundance (**Figure 5A**). Partial redundancy analysis showed that the bacterial community, MGEs and environmental factors separately explained 19.23, 19.87, and 9.4% of the variation of ARG profiles, respectively. Interaction between

all variables explained 4.63% of the variation, and 42.49% of the variation was unexplained.

### Correlations Between ARGs and MGEs

The absolute abundance of ARGs was positively correlated with MGEs (**Supplementary Figure S8**), and the normalized abundance of MGEs was significantly correlated with the richness and normalized abundance of ARGs (**Supplementary Table S4**). The network analysis further indicated that aminoglycoside, beta-lactamase, MLSB, and tetracycline resistance genes were significantly correlated with MGEs (**Figures 3C,D**), and the co-occurrence pattern revealed that Tp614 and tnpA-03 play the most important roles in the network (**Supplementary Figure S7B**). However, the integrase gene intI-1 only showed a significant correlation with multidrug and vancomycin resistance genes (**Supplementary Table S5**), but the absolute abundance of transpose gene tnpA-05 showed a strong correlation to all ARG types except vancomycin resistance genes (**Supplementary Table S6**). Additionally, the transposase genes IS613, tnpA-03, tnpA-04, and tnpA-07 also showed a significant correlation with some types of ARGs. Partial redundancy analysis showed that 19.87% of the variation of ARG profiles was explained by MGEs (**Figure 5B**), indicating that MGEs were an important contributor to the change

size times immigration, R 2 indicates the fit to the neutral model.

of ARG profiles during our study period. These results suggested that MGEs played a crucial role in the change, accumulation and dissemination of ARGs in this urban landscape reservoir.

# DISCUSSION

### The Contrasting Temporal Patterns Between Bacterial OTUs and ARGs

Our high-frequency time series-based study revealed clear differences in the temporal patterns of microbial communities and ARG profiles, along with the decoupling of the temporal dynamics of bacterial taxonomy and antibiotic resistance function in an urban waterbody. The composition of ARG subtypes varied frequently over the year but was relatively stable at the ARG type level. None of the ARG subtypes was present in all 51 samples, and no distinct seasonal patterns in ARG profiles could be determined, but microbial taxonomic community had significant seasonal pattern.

ARGs in this urban reservoir were detected at relatively high abundance and richness. The diversity of ARGs in this study was greater than that reported from Chinese suburban/rural lakes/reservoirs (Liu et al., 2018). The high abundance and richness in this waterbody may be driven by the combined human activities of sewage discharge and agricultural/urban runoff, as it is generally known that these activities increase the number of ARGs in waterbodies (Zhu et al., 2017; Ahmed et al., 2018; Liu et al., 2018). Moreover, aminoglycoside and multidrug resistance genes dominated the abundance of ARGs in this study, while aminoglycoside and beta-lactamase resistance genes dominated the types or richness of ARGs. In contrast, Jia et al. (2017) found in swine wastewater and downstream water that tetracycline and aminoglycoside resistance genes were the dominant ARGs, perhaps due to the different antibiotic selection pressures in different environments. However, the primer limitation of this study may result in the bias of detected ARGs, as the PCRbased method only covered limited number of ARG types. Additionally, a majority of ARGs investigations have been carried out in municipal wastewater treatment plants. An et al. (2018) showed significant changes of ARG profiles based on samples just collected in February and August, while another study based on monthly samples for 1 year did not find obvious seasonal variation in methicillin-resistant gene although it varied between

months and sampling sites (Börjesson et al., 2009). Yang et al. (2013) suggested that some types of ARGs in activated sludge showed seasonal fluctuations during a 4-year investigation, while most types of ARGs did not show seasonal fluctuations, but these results were based on samples collected only in summer and winter. However, according to our study, ARG subtypes were extremely dynamic and some of them frequently showed rises and falls over short-time periods (i.e., in just one week) but did not show seasonal pattern, so the previous low-frequency based studies may not completely and effectively reflect actual seasonal dynamic patterns of ARGs at fine temporal scale. Overall, our study based on high-frequency sampling over long timescales provide better understanding of the detailed temporal pattern of ARG profiles.

Our data revealed that some bacterial OTUs (e.g., Planktothrix) fluctuated rapidly, while major phyla remained stable at a broad taxonomic level. The microbial community showed an obviously seasonal variation and a trend of resilience during 1 year (12 months). In the microbial community, Proteobacteria and Actinobacteria were the dominant phyla, consistent with previous studies in drinking water system (Jia et al., 2015) and rivers (Jia et al., 2017). Additionally, the high relative percentage of Cyanobacteria was previously reported to have significant impacts on the variation of ARGs (Guo et al., 2018). Wang et al. (2018) showed that Proteobacteria, Actinobacteria and Cyanobacteria in the environment can carry some ARGs, confirming the high normalized abundance of ARGs and the vital roles of these bacteria in the dissemination of ARGs. Furthermore, we found a significant seasonality in the microbial taxonomic community, and the Procrustes test indicated the asynchronous temporal patterns between microbial OTUs and ARGs compositions (**Figure 2**). Previous studies have shown that ARGs strongly correlate with microbial communities in antibiotic-rich environments including human gut (Feng et al., 2018) and landfill leachate (Zhao et al., 2018). Recently, Ma et al. (2017) indicated that ARGs and bacterial communities did not show a strong correlation in drinking water, which is similar to the result of our study under low or medium antibiotic pressure. Altogether, the weak correlation between ARGs and bacterial communities in this study may be due to the relatively low antibiotic and heavy metal selection pressure in the natural environment.

# The Co-occurrence Pattern of Bacterial OTUs and ARGs

Previous studies have confirmed that network analysis is a useful and reliable method to decipher the co-occurrence patterns between bacterial taxa and ARG subtypes, and to track potential hosts of ARGs (Barberán et al., 2012; Li et al., 2015; Ma et al., 2017; Feng et al., 2018). The module hubs in the network function as keystone of the module (Faust and Raes, 2012). All

module hubs in our co-occurrence network were affiliated to ARGs, implying that ARGs were important in maintaining the structure of the network and that HGT of ARGs was frequent among bacteria in the microbial community. The high frequency of aminoglycoside resistance genes in the network may be due to their high abundance and richness in the community. However, the abundance and richness of tetracycline resistance genes were relatively low in our studied system, but they connected tightly with various bacterial taxa. This suggests that the tetracycline resistance genes were shared among multiple bacterial taxa (Li et al., 2015). An increasing number of studies have reported the widespread distribution of tetracycline resistance genes (Czekalski et al., 2014; Zhou et al., 2017), as this antibiotic has been widely used in China (Hvistendahl, 2012).

In the co-occurrence network, Proteobacteria, Bacteroidetes and Firmicutes were the most prevalent phyla. These taxa were predicted as possible ARG hosts, and/or they originated from a similar source including upstream wastewater and surface runoff, suggesting the potential dissemination risk of ARGs with these taxa. In fact, Proteobacteria include a wide variety of pathogens and occur frequently as opportunistic pathogens, and have been found in soils that contained almost all major ARG types (Forsberg et al., 2014). Moreover, Jang et al. (2018) found that Bacteroidetes can be potential hosts of the intI-1 gene. In this study, Bacteroidetes were significantly correlated to all major ARG types except sulfonamide resistance genes (**Figure 3** and **Supplementary Table S2**), suggesting the potential transmission of ARGs to human pathogenic bacteria. The Firmicutes were strongly correlated to aminoglycoside resistance genes in sewage sludge (Su et al., 2015), and the ARGs were shown to frequently switch hosts from Firmicutes to other bacterial phyla (Wang et al., 2017). In our study, Firmicutes were not the dominant bacteria but may carry diverse ARGs and further cause ARGs pollution to different environments or ecosystems. Nevertheless, the co-occurrence patterns revealed by this study were based on correlation, so the exact ARG carrying host and transmission pathway of ARGs need to be verified by further study.

### Assembly Mechanisms Driving the Dynamics of ARGs

Our results indicated that the variation in ARGs could be explained by the bacterial community, MGEs and environmental factors. Environmental factors including salinity, total carbon, nitrate nitrogen, total phosphorus and precipitation influenced the temporal variation of ARGs, while precipitation was the strongest factor. A recent study reported that nutrients explained a large part of ARGs variation and that the nutrients combined with other factors could drive the distribution of ARGs (Zhao et al., 2017). The seasonal stormwater runoff induced by rainfall can increase the absolute abundance of ARGs and MGEs in an aquatic environment (Chen et al., 2019) and could influence the ARGs load through increased input from runoff and soil resuspension (Di Cesare et al., 2017). In accordance with previous studies, we found that precipitation displayed a significant influence on the ARG composition and structure, and was positively correlated with the richness and normalized abundance of ARGs, supporting the view that rainfall contributed to the loading of ARGs in the river and reservoir waters in urban region. In addition, rainfall significantly influenced pathogenic bacterial abundance because stormwater was able to significantly contribute to the occurrence and elevated concentration of pathogenic bacteria in subtropical waters (Ahmed et al., 2018). Our results further showed that turbidity was strongly correlated to richness and diversity of bacteria and ARGs. A turbidity increase is normally associated with rainfall-induced contamination, suspended particulate matter and sediment disturbance. This may further drive changes in bacterial community and ARG profiles. It has been demonstrated that turbidity was significantly correlated to bacterial diversity, reiterating the influence of runoff on the bacterial community (Peter and Sommaruga, 2016). Taken together, our results indicated that precipitation had both direct and indirect effects on bacterial taxonomic communities and resistomes, supporting the idea that stormwater with agriculture/urban runoff, with contaminants from the urban region could increase loadings of bacterial ARGs thereby increasing the risk of ARGs dissemination (Garner et al., 2017). In contrast, we did not find any significant relationship between water temperature and the ARGs composition or abundance, corroborating the fact that the ARGs in this subtropical urban reservoir did not have a significant seasonal pattern. This is inconsistent with a previous finding that water temperature was a potential factor driving the dynamics of ARG profiles (Zheng et al., 2018), and this inconsistency may be attributed to a difference in sampling frequency and climate conditions, as the range of temperature was narrower in this subtropical reservoir compared to other studies in temperate zone. In the Zheng et al. (2018) study, water samples were only collected in March, June, September, and December, and their river showed a wider temperature range than our studied reservoir. In our study, the 42.49% of ARGs variance remained unexplained, may partly be due to the environmental factors that were not considered, such as anthropogenic factors and the co-selection of antibiotics and heavy metals (Stepanauskas et al., 2006; Liu et al., 2018; Chen et al., 2019).

The bacterial community has been identified as one of the key drivers that shape the ARG profiles in antibioticrich environments (Forsberg et al., 2014; Su et al., 2015). However, our results of pRDA, Procrustes test and Spearman's relationship all showed that microbial community was weakly correlated with ARG profiles. More interestingly, our results indicated that the neutral process played a more important role in the assembly of microbial taxonomic communities than in ARG profiles, as the neutral model only explained a minor part (7.9%) of the variation in ARG profiles. Therefore, dispersal and stochastic processes appeared to have more influence on bacterial taxonomic communities compared to ARG profiles. This may be one of the reasons why microbial OTUs and ARGs communities exhibited a decoupling change over 1 year. These results indicated that the vertical transmission

(mother-to-child transmission) of ARGs in this environment may not play a leading role.

Our study showed that MGEs could partly explain the variation of ARG profiles in the urban reservoir, and the absolute abundance of MGEs showed a significant and positive correlation with ARGs. Further, MGEs that were always abundant in our samples, could potentially enhance the probability of HGT of ARGs. Previous studies reported that environmental contaminants at low concentration have a potential role in the dissemination of ARGs by promoting HGT, such as antibiotic and non-antibiotic pharmaceuticals (Andersson and Hughes, 2014; Wang et al., 2019). In this study, the rainfall combined with possible contaminants and MGEs may largely facilitate the proliferation and spread of ARGs among different bacterial taxa. These results suggest that MGEs played an important role in the temporal variation of ARGs in the urban waterbody, implying the potential risk that the ARGs may be transferred by MGEs from aquatic environments to human pathogens (Martínez et al., 2015).

### CONCLUSION

Our weekly study revealed the asynchronous succession patterns of bacterial taxonomic communities and ARG profiles, indicating a decoupling dynamic between bacterial taxonomy and function on a temporal scale in a subtropical urban reservoir. The bacterial taxonomic community was significantly correlated with water temperature and exhibited seasonality, while ARG profiles was strongly correlated with precipitation and did not show distinct seasonality. Precipitation had both direct and indirect effects on the temporal patterns of bacterial community and ARG profiles, reflecting the effect of stormwater runoff on temporal dynamics of ARG profiles. The diversity and composition of bacterial taxonomic communities appeared to be more responsive to environmental variables than ARGs, and the composition of ARG profiles appeared to be governed less by stochastic processes. Additionally, MGEs were positively correlated with ARGs and had significant contribution to the variation of ARGs community, suggesting the potential risk of HGT. This study may help to better understand the temporal dynamic of ARG profiles in urban water environments at fine temporal scale, and highlight that further high-frequency and long-term time series research should be undertaken to

### REFERENCES

reveal the specific drivers and deep mechanisms of ARGs assembly and dynamics.

### DATA AVAILABILITY

The raw data for 16S rRNA gene analysis in this study have been deposited in the NCBI SRA database under the BioProject number PRJNA510463 and the accession number SRP173857.

### AUTHOR CONTRIBUTIONS

JY designed the experiments. PF, FP, and XG carried out the sample collection and determined the environmental parameters. PF, FP, and PX conducted the PCR, highthroughput qPCR, high-throughput sequencing, and bioinformatics. PF and JY wrote the first draft of the manuscript. All authors contributed to and have approved the final manuscript.

### FUNDING

This work was funded by the Strategic Priority Research Program of the Chinese Academy of Sciences (XDA23040302), the National Natural Science Foundation of China (91851104), the Key Laboratory of Urban Environment and Health, Chinese Academy of Sciences (KLUEH-C-201801), and the Xiamen Municipal Bureau of Science and Technology (3502Z20172024 and 3502Z20171003).

### ACKNOWLEDGMENTS

We thank Prof. David M. Wilkinson for his comments on the manuscript.

### SUPPLEMENTARY MATERIAL

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

river-reservoir system. Water Res. 149, 179–189. doi: 10.1016/j.watres.2018. 11.007

distribution systems. Environ. Pollut. 213, 119–126. doi: 10.1016/j.envpol.2016. 02.013

**Conflict of Interest Statement:** 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.

Copyright © 2019 Fang, Peng, Gao, Xiao and Yang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# The Relative Abundance of Benthic Bacterial Phyla Along a Water-Depth Gradient in a Plateau Lake: Physical, Chemical, and Biotic Drivers

*Kaiyuan Wu1,2† , Wenqian Zhao2,3† , Qian Wang2 , Xiangdong Yang2 , Lifeng Zhu3 , Ji Shen2 , Xiaoying Cheng1 \* and Jianjun Wang2,4 \**

### *Edited by:*

*Haihan Zhang, Xi'an University of Architecture and Technology, China*

### *Reviewed by:*

*Jie Wang, China Agricultural University, China Man Kit Cheung, The Chinese University of Hong Kong, China*

### *\*Correspondence:*

*Xiaoying Cheng<EMAIL_ADDRESS>Jianjun Wang<EMAIL_ADDRESS>
*† These authors have contributed equally to this work*

### *Specialty section:*

*This article was submitted to Aquatic Microbiology, a section of the journal Frontiers in Microbiology*

*Received: 30 April 2019 Accepted: 18 June 2019 Published: 10 July 2019*

### *Citation:*

*Wu K, Zhao W, Wang Q, Yang X, Zhu L, Shen J, Cheng X and Wang J (2019) The Relative Abundance of Benthic Bacterial Phyla Along a Water-Depth Gradient in a Plateau Lake: Physical, Chemical, and Biotic Drivers. Front. Microbiol. 10:1521. doi: 10.3389/fmicb.2019.01521*

Frontiers in Microbiology | www.frontiersin.org 1 July 2019 | Volume 10 | Article 1521

*1 School of Environment and Civil Engineering, Jiangnan University, Wuxi, China, 2 State Key Laboratory of Lake Science and Environment, Nanjing Institute of Geography and Limnology, Chinese Academy of Sciences, Nanjing, China, 3 School of Biological Sciences, Nanjing Normal University, Nanjing, China, 4 University of Chinese Academy of Sciences, Beijing, China*

Water-depth biodiversity gradient, one of the typical biogeographical patterns on Earth, is understudied for bacteria in freshwater ecosystems, and thus left the underlying mechanisms poorly understood especially for benthic bacteria. Here, we investigated the water-depth distribution of surface sediment bacterial phyla and their driving factors in Lake Lugu, a plateau lake in Southwest China. Our results revealed that the relative abundance of 11 dominant bacterial phyla showed various water-depth patterns, such as increasing, decreasing, hump-shaped, and U-shaped patterns. These patterns across phyla were consistent with their different niche positions of water depth, while the occupancy-abundance relationships were not dependent on phylum attributes. Consistently, phylum abundance was best explained by water depth; other physical and chemical factors, such as metal ion concentrations, SiO2, and pH, can also explain the variations in some bacterial phyla. Chemical variables were the main drivers of the dominant bacterial phyla. However, biotic variables also showed substantial importance for some phyla, such as Planctomycetes, Actinobacteria, and WS3. This work could provide new insights into the general water-depth patterns and underlying mechanisms of the relative abundance of bacterial phyla in freshwater ecosystems.

Keywords: bacteria, phyla, relative abundance, water depth, biotic effects, abiotic effects

# INTRODUCTION

Biogeographical patterns of biodiversity along gradients, such as those of latitude, elevation, and water depth, are among the most widely studied topics in ecology (Zintzen et al., 2017), and researchers have described the latitudinal and elevational patterns in biodiversity in a variety of ecosystems across the globe (Fuhrman et al., 2008; Wang et al., 2011; Shi et al., 2014). Similar to the more commonly studied latitudinal and elevational gradients, water-depth gradient show complex changes in various environmental attributes and thus impose strong environmental filtering on the aquatic microbial community (Bryant et al., 2012). The phylogenetic and physiological diversity of microbial phyla, especially those of bacteria, is considerably greater than that of animal and plant phyla (Prosser et al., 2007; Yeh et al., 2019). Furthermore, benthic bacteria play an important role in organic matter transformation and in the biogeochemical cycling of major elements such as nitrogen, phosphorus, sulfur, and iron (Martins et al., 2011; Cheng et al., 2014). Knowledge of benthic bacterial phyla distribution in aquatic ecosystems is of great significance and can provide novel perspectives regarding water-depth biodiversity patterns among microorganisms and their potential drivers and further contribute to the sustainable management of water resources. Nonetheless, as most studies focus on marine waters (Agogue et al., 2011; Bryant et al., 2012; Walsh et al., 2016), the general patterns and driving mechanisms of bacterial communities along water-depth gradients have been less studied for freshwater ecosystems.

In recent years, the niche position and niche breadth hypotheses have been used to investigate how niche characteristics account for species distribution and abundance (Brandle and Brandl, 2001). Among these two referenced niche parameters, niche position indicates the availability of habitat and resources (Teittinen et al., 2018). Some species have marginal niche positions at the very edges of the environmental ranges, while other species have nonmarginal niche positions that exist under average environmental conditions (Vilmi et al., 2019). Niche position, rather than niche breadth, is the chief predictive factor of species distributions, such as those of stream insects (Heino and de Mendoza, 2016). In addition, previous research has provided strong evidence that the occupancy and abundance relationships of species tend to be positively and often strongly linked with each other in a variety of ecosystems (Gaston et al., 2000; Foggo et al., 2007; Zuckerberg et al., 2009). Niche position could be a strong predictor of species variation in terms of occupancy and abundance. For instance, niche position is negatively associated with mean abundance in diatom and macroinvertebrate species, and the occupancy-abundance relationship clearly differs across taxonomic groups (Vilmi et al., 2019). At the community level, niche position and occupancy are related to the abundance distribution, but whether such a relationship exists at phylum level remains unexplored.

Traditional niche-based theory asserts that the relative abundances of species are determined by abiotic and biotic factors, such as environmental conditions, habitat heterogeneity, and species interactions (Dumbrell et al., 2010; Liao et al., 2016). Most studies have shown that local environmental conditions, such as physical and chemical attributes, can affect microbial community composition and species diversity in lake ecosystems and the biogeochemical processes they mediate (Staley et al., 2015; Qin et al., 2016). The bacterial community can be shaped by a variety of environmental factors, such as water depth, nitrogen, phosphorus, pH, and pollution (Brown et al., 2009; Haller et al., 2011; Song et al., 2012; Liu et al., 2015b; Zhang et al., 2019). In addition, biotic interactions could also affect microbial community structure. However, rarely considered are the important ecological processes associated with biotic variables, such as the community composition and diversity of other organisms that feed on bacteria (Langenheder et al., 2017). Ecosystems are constructed around interaction webs, such as predation, mutualism, competition, and hostparasite interactions, that connect every species to many others at an array of spatial scales (Estes et al., 2011; Amin et al., 2012). For example, viral abundance is typically strongly correlated with bacterial abundance in marine systems (Clasen et al., 2008; Hanson et al., 2017). In the soil, fungi and bacteria have an effect on the formation of each other's community structure (Singh et al., 2009). Moreover, biotic attributes such as the composition of other taxonomic groups are generally better predictors of species diversity in streams than other environmental variables (Johnson and Hering, 2010).

Here, we explored the water-depth patterns and underlying drivers of the relative abundance of bacterial phyla in the surface sediments of Lake Lugu, a typical deep plateau lake system located in Southwest China. This lake has a water depth of 93.5 m and an area of 50.5 km2 . We had the four following objectives. First, we examined the water-depth patterns in terms of the relative abundance of benthic bacterial phyla. Second, we expected that niche positions at the phylum and species levels are related and also the relationship between the occupancy and abundance of each phylum. Finally, we explored the main factors driving the above water-depth patterns in bacterial phyla and further quantified the relative importance of physical, chemical, and biological drivers.

### MATERIALS AND METHODS

### Study Region

The studied Lake Lugu, with an elevation of 2,685 m, is located in the Yunnan Province, Southwest China (27°41′– 27°45′N, 100°45′–100°50′E). It is one of the deepest freshwater plateau lakes in the region. This lake has spatial characteristics that include high connectivity and a small area, allowing species to be freely distributed at different locations with no apparent barriers to dispersal. Macrophytes were partly absent from the lake, allowing for the collection of surface sediments in shallow waters. In August 2010, we collected 37 surface sediment (~0–1 cm) samples along the water-depth gradient from 0 to 93.5 m. At each site, three sediment cores with a 6-cm diameter were retrieved for surface sediments. The sample collection process was described in detail in a previous study (Wang et al., 2012b).

### Bacterial Community Analysis

Bacterial analyses were performed according to the previous literature (Wang et al., 2012a). Briefly, genomic DNA was extracted from surface sediment samples using the phenol chloroform method (Zhou et al., 1996). Bacterial tag-encoded FLX amplicon pyrosequencing was conducted as described previously (Dowd et al., 2008). Bacterial 16S rRNA genes were amplified using the 27F primer (5′ GAG TTT GAT CNT GGC TCA G 3′) with the 454 Life Sciences A sequencing adapter, and the modified 519R primer (5′ GTN TTA CNG CGG CKG CTG 3′) with a 8-bp barcode sequence and the 454 Life Sciences B sequencing adapter. The PCR amplification of a single 35 cycle was performed in a ABI9700 thermocycler (ABI, Foster City, USA) using the program 95°C for 2 min; 25 cycles of 95°C for 30 s, 55°C for 30 s, 72°C for 30 s; 72°C for 5 min; finally kept at 10°C. Triplicate positive PCR products were pooled and purified with AxyPrep DNA Gel Extraction Kit (Axygen, USA). Amplicon sequencing was performed based upon the manufacturer protocols (Roche Applied Science, Indianapolis, IN) for Titanium sequencing on FLX-titanium platform (Dowd et al., 2008). The purified amplicons were pooled at equal molality, and then sequenced using a Roche 454 FLX pyrosequencer (Roche, Switzerland).

All 16S rRNA pyrosequencing reads were analyzed using QIIME, version 1.9.0 (Caporaso et al., 2010b). We imported barcoded 16S rRNA gene sequences and removed the primers, demultiplexed reads, and then filtered them according to Phred quality scores. Quality criteria were a minimum sequence length of 200 bp, a maximum sequence length of 900 bp, and a minimum average quality score of 25. With such quality control, we allowed few ambiguous bases or mismatches in the primer sequence and no barcode errors, and a maximum homopolymer length of 6 bp.

Sequences were denoised with the Denoiser algorithm (Reeder and Knight, 2010), clustered into operational taxonomic units (OTUs) at 97% similarity level with the seed-based UCLUST algorithm (Edgar, 2010). We used the longest sequence in a cluster as the representative sequence for that OTU. Singletons, i.e., OTUs with only one read in the entire dataset, were removed and chimeras were excluded using ChimeraSlayer (Haas et al., 2011). After chimeras were removed *via* ChimeraSlayer, representative sequences from each OTU were aligned to the Greengenes v13.8 imputed core reference alignment (DeSantis et al., 2006) using PyNAST (Caporaso et al., 2010a). The taxonomic identity of each representative sequence was determined using the RDP Classifier (Wang et al., 2007) and chloroplast and archaeal sequences were removed. Before the subsequent analysis, bacterial community were rarefied at 1,139 sequences to avoid the variation in abundance or sampling intensity biased. The generated sequences can be found in figshare (https://doi.org/10.6084/ m9.figshare.8052788).

### Abiotic and Biotic Variables

We evaluated several environmental characteristics important to the biological community in the surface water of 0.5 m and bottom water near the sediment-water interface at each site (**Supplementary Table S1**). We analyzed the temperature, pH, conductivity, total nitrogen, total phosphorus, dissolved oxygen (DO), HCO3 - concentration (HCO3·water), and silica of the surface water, and measured the temperature, pH, conductivity, and DO of the bottom water.  The measured metal ions included Al, Ba, Be, Ca, Co, Cr, Cu, Fe, Li, K, Mg, Mn, Na, Ni, Pb, Sr., Ti, V, and Zn. The grain size was divided into five classes: <4, 4–16, 16–32, 32–64, and >64 μm. A detailed description of the methods for the measurement and calculation of the abiotic variables was described in a previous study (Wang et al., 2012b).

In terms of the biotic variables, diatoms and chironomids were classified and identified to the lowest possible taxonomic level (typically to the species level) as shown in previous studies (Wang et al., 2012b; Zhang et al., 2013). We used the following biotic factors as predictor variables: (1) the chlorophyll *a* concentration of the surface and bottom water, (2) the biomass of diatoms and chironomids, and (3) the richness of diatoms and chironomids. The concentration of chlorophyll *a* was the representative of the phytoplankton biomass. The biomass of diatoms and chironomids was indicated by the number of diatoms per gram of wet sediment and the count of head capsules per gram of wet sediment, respectively.

### Statistical Analyses

First, we only included the highly abundant bacterial phyla by considering species occurrence percentages larger than 85%. We explored the relationships between water depth and the relative abundance of bacterial phyla with linear or quadratic models. The most appropriate model was selected based on the lowest value of Akaike's information criterion (Yamaoka et al., 1978).

Second, we calculated the niche position of the water depth for both the species and phylum levels. For each species or phylum, we estimated the weighted water-depth position by averaging the product of their water-depth values and their abundances across all samples (Wang et al., 2013). We selected the linear or quadratic model based on the lower value of Akaike's information criterion and explored the relationship of niche position between the phylum and species levels. In practice, the proportion of sites was considered as "occupancy," and the mean abundance at occupied sites was considered as "abundance." The relationship between the occupancy and abundance of bacterial phyla was examined by linear or quadric models, and the better model was selected based on the lower value of Akaike's information criterion (Yamaoka et al., 1978).

Third, to identify the important factors affecting each phylum, we applied random forest model (Feld et al., 2016) and redundancy analysis (RDA). Statistical dependence between the explanatory variables was assessed using Pearson's correlation coefficients, and variables with high correlation coefficients (Pearson *r* > 0.7) were excluded from the models. By performing principal component analysis, the 19 metal ion concentrations were reduced to the first two axes of the principal components analysis, which thus represented the environmental variables that reflected the geochemical factors. This process was performed to decrease the degrees of freedom so that they were lower than the number of sampled sites. The other measured variables were used as environmental variables without principal component analysis. For the random forest analysis, we used techniques such as cross-validation to prune the 2,000 trees to an optimal size (Prasad et al., 2006). The importance of a predictor variable was determined by its frequency of selection weighted by a measure of improvement of the model given each split and averaged across all trees (contributions were scaled to sum to 100).

Furthermore, RDA was used to examine the potential explanatory variables of community composition at the phylum level. We used RDA with Hellinger-transformed abundance data of the bacterial community, as this transformation makes complex data with numerous zero values more suitable for the analysis with linear methods (Legendre and Gallagher, 2001).

Finally, we applied variation partitioning to quantify the relative importance of the main environmental drivers using linear model (Anderson and Cribble, 1998; Legendre, 2008). We categorized the explanatory variables into three groups: physical, chemical, and biological components. The physical component included water depth, conductivity of surface water, porosity, and grain size. The chemical component included surface water pH, total phosphorus, total nitrogen, SiO2 and HCO3, bottom water DO and pH, and surface sediment LOI, grain size, PC1 and PC2 of metals. For biological component, we considered the richness and biomass of diatoms and chironomids, and chlorophyll *a* of surface and bottom water. We generated three explanatory matrix models and estimated the proportions of variation in the relative abundance of the bacterial phyla explained by these three components. For each component, significant explanatory variables were selected by forward selection against the Hellinger-transformed abundance phylum data with 9,999 permutations.

We did not apply more statistical analyses because it is a challenge to differentiate the direct or indirect effects of underlying drivers on the relative abundance of bacterial phyla based on field observations. These above analyses were conducted with vegan V2.5-5 (Oksanen et al. 2013) and randomForestSRC V2.9.0 (Ishwaran and Kogalur 2014) in the R environment.

### RESULTS

In total, we selected 11 major bacterial phyla based on their relative abundance, including Proteobacteria, Nitrospirae, Chlorobi, Chloroflexi, and Bacteroidetes. Among these phyla, Proteobacteria was the most dominant, followed by Bacteroidetes, Chloroflexi, and Nitrospirae (**Table 1**). The relative abundance of most phyla, seven out of 11, exhibited significant (*p* < 0.05) and clear water-depth patterns (**Figure 1**), such as increasing, decreasing, hump-shaped, and U-shaped trends. For example, the relative abundance of Bacteroidetes decreased with water depth, while that of Nitrospirae increased toward deep water (**Figure 1**). In addition, Acidobacteria, Planctomycetes, and

TABLE 1 | Abbreviations for the abundant phyla.

*The Table includes the most abundant bacterial phyla that together contributed up to 85% of the total sequence read. Relative abundance represents the average value of each phylum across all sites.*

modeled with linear and quadratic models. The better model was selected based on the lower value of Akaike's information criterion. For better visualization, species relative abundance was scaled as mean = 0 and SD = 1.

Gemmatimonadetes showed hump-shaped trends along the water-depth gradient (**Figure 1**). Interestingly, the phototrophic phyla Chlorobi and Chloroflexi showed U-shaped patterns.

At the phylum level, the niche position of water depth showed great variations, ranging from 26 to 48 m (**Figure 2A**). The phylum-level niche positions had significant (*p* < 0.05) positive correlations with those at the species level (**Figure 2A**). For the whole bacterial community, there was a significant relationship between occupancy and abundance (**Figure 2B**). This is consistent with the finding at the phylum level, which also showed rather strong and positive occupancy-abundance relationships (**Figure 2B**). However, there were no clear differences among phyla in terms of the relationship of occupancy and mean abundance.

The random forest analyses showed that water depth was the strongest variable correlating with the relative abundance of the most bacterial phyla. Other environmental variables, such as metal ion concentrations, SiO2, pH, total phosphorus, DO, and conductivity also have important effects on bacterial phyla (**Figure 3A**). Interestingly, biological factors such as chironomid richness and chlorophyll *a* concentration were also important for some bacterial phyla (**Figure 3A**). Such findings were confirmed by RDA, which shows that water depth was the strongest factor (*p* < 0.001) related to the distribution of bacterial phyla, while biological attributes, such as diatoms biomass, also showed significant (*p* = 0.038) correlations (**Figure 3B**).

In the variation partitioning analyses with the three variable components, the pure effect of chemical component was the higher for bacterial phyla rather than the pure effects of physical and biological components. The relative abundance of most phyla, such as Planctomycetes (29.9%), WS3 (22.1%), and Verrucomicrobia (21.0%), was well explained by the pure

FIGURE 2 | The relationship of niche position between the phylum and species levels (A), and the relationships between the occupancy and abundance of each phylum (B). The shapes in (A) represent various distributional patterns: ■, hump-shaped pattern; ●, decreasing pattern; ▲, increasing pattern; ♦, U-shaped pattern. (B) Relationships between occupancy and mean local abundance for bacterial phyla. The points in (B) represent the different species. Linear or quadratic models were selected according to a lower value of Akaike's information criterion.

FIGURE 3 | The abiotic and biotic factors related to the relative abundance of bacterial phyla. These factors were identified with random forest (A) and redundancy analysis (RDA, B). For abiotic factors, we considered water depth (depth), total phosphorus (TP.water), total nitrogen (TN.water), conductivity (Con.sur) and concentrations of - HCO3 (HCO3·water) and SiO2 (SiO2·water), pH (pH.bot and pH.sur), dissolved oxygen (DO.bot), sediment porosity, loss-on-ignition (LOI), the first two axes of metal ion principal component analysis (Comp.1 and Comp.2) and grain size. Grain size: <4 μm (GSL4), 4–16 μm (GS4\_16), 32–64 μm (GS32\_64), and > 64 μm (GSL64). Biotic factors included the biomass (Dia.bio and Chi.bio) and richness (Dia.rich and Chi,rich) of diatoms and chironomids, and chlorophyll *a* (Chla.bot and Chla.sur). For RDA, the abiotic variables were automatically selected based on Monte Carlo permutation tests (999 permutations).

effects of the three components in predicting the relative abundance are shown, but not the joint effects or unexplained variances. An alternative version of this

effect of chemical component. This is especially true for Verrucomicrobia and Gemmatimonadetes, which were only explained by the pure effect of chemical component (**Figure 4**). Physical and biological variables had effects on seven and six bacterial phyla, respectively. For Nitrospirae, physical variables were especially important, which explained 16.3% of its variation. The pure effect of biotic variables accounted for 8.1% variations of the relative abundance of Planctomycetes. For Proteobacteria, only 0.1% was explained by biological variables. Four bacterial phyla can be explained by the combination of pure effects of biological, chemical, and physical variables.

figure showing the unique and shared variance of each group can be found in Supplementary Figure S1.

# DISCUSSION

A long-standing task in ecology has been to explain biological distribution patterns and the drivers underlying these patterns. In recent years, the water-depth patterns of microbial community distribution have been extensively investigated in aquatic systems. However, most studies have focused on marine environments (Smith and Brown, 2002; Wang et al., 2013), while few studies have been conducted in freshwater ecosystems, such as lake sediments (Haglund et al., 2003; De Wever et al., 2005; Yu et al., 2014). The main aim of our study was to explore the distributional patterns of bacterial abundance at the phylum level and to analyze the abiotic and biotic variables explaining these patterns. Our results highlight four main findings. First, the bacterial phyla exhibited various patterns along the water-depth gradient, some of which are rarely reported, such as the hump-shaped and U-shaped patterns. Second, the occurrence of these patterns across phyla might be explained by the variation in their niche positions, while the occupancy-abundance relationships were not dependent on phylum attributes. Third, among all measured exploratory variables, water depth was the most important predictor of the relative abundance of bacterial phyla. Finally, although the environmental variables were shown to have pivotal roles in terms of the relative abundance of bacterial phyla, biotic variables could also substantially explain the abundance of some phyla, such as Planctomycetes, Actinobacteria, and WS3.

Proteobacteria was the most dominant phyla in the sediments of this lake, followed by Bacteroidetes, Chloroflexi, and Nitrospirae. This is consistent with previous studies that Proteobacteria and Bacteroidetes are the most abundant phyla in sediments (Diao et al., 2017). These phyla have frequently been shown to dominate freshwater bacterial communities (Dai et al., 2016). At present, investigations of the spatial pattern of bacterial relative abundance in freshwater lake sediments indicate that bacterial abundance shows clear spatial variations (Bouzat et al., 2013; Liu et al., 2015a). However, the distributional patterns of bacterial abundance along water-depth gradients, especially among bacterial phyla, have been less explored. Bacterial abundance generally decreases with water depth in freshwater systems (Liu et al., 2016). Our results reveal various water-depth patterns in terms of the relative abundance of bacterial phyla, among which some patterns have been rarely reported so far. For instance, there were significantly (*p* < 0.05) decreasing patterns for Bacteroidetes, increasing patterns for Nitrospirae, and hump-shaped patterns for Acidobacteria, Planctomycetes, and Gemmatimonadetes. The decreasing and increasing water-depth patterns of the relative abundance of bacterial phyla are consistent with the patterns found in earlier studies along other environmental gradients, such as elevations. For instance, the relative abundance of Actinobacteria, Bacteroidetes, and Deltaproteobacteria increases with elevation in stream biofilms, while that of Alphaproteobacteria decreases with elevation (Wang et al., 2012a). However, the U-shaped patterns of Chlorobi and Chloroflexi were unexpected because such patterns are extremely rarely observed in nature, and these two phyla have phototrophic capability (Gupta, 2010). The two phyla are commonly present in the surface waters and bottom waters of deep lakes and oceans (Newton et al., 2011). The phylum Chlorobi is resistant to very low light conditions (Freed et al., 2019) and primarily exists in the anoxic aquatic environment of stratified lakes (Gupta, 2010). Chloroflexi, however, contains both photosynthetic and nonphotosynthetic members (Bryant and Frigaard, 2006), but metabolic insights suggest that a primarily heterotrophic lifestyle boosted by light-driven energy generation *via* aerobic anoxygenic phototrophy occurs in some groups (Mehrshad et al., 2018). It should be noted that alternative explanations may be needed for such U-shaped patterns of Chlorobi and Chloroflexi.

These water-depth patterns may be explained by the niche positions of the bacterial phyla. Consistent with the various water-depth patterns of relative abundance, the niche position also obviously differed among phyla. For instance, the phyla with decreasing water-depth patterns had a lower niche position, while those with increasing and hump-shaped patterns had a higher niche position. In other words, niche position may have a certain impact on bacterial phylum distributional patterns. This phenomenon has also been reported for other taxonomic groups, such as the distribution of stream insects, the predictability of which depends on niche position (Heino and de Mendoza, 2016). Furthermore, niche position could be an important predictor of the regional occupancy and local abundance of diatoms and insects, such as those in streams (Rocha et al., 2018). However, the occupancy-abundance relationship showed a relatively similar trend for all phyla. This result may have occurred because of the taxonomic relatedness among bacteria, as taxonomic relatedness has only minor influence on occupancy and abundance, as observed among the diatoms and macroinvertebrates across a set of lakes (Heino and Tolonen, 2018). Hence, the occupancy-abundance relationship could not explain the variation in water-depth patterns across phyla.

Physicochemical variables were important in explaining the relative abundance of bacterial phyla along the water-depth gradient. Among them, water depth was the most important environmental factor affecting the relative abundance of bacterial phyla in this lake ecosystem, which is consistent with previous studies showing that water depth has important effects on bacterial community composition (Yokokawa et al., 2010; Zhang et al., 2015). However, water depth could also be associated with many other physiochemical variables, including temperature, pressure, productivity, and nutrient availability (Smith and Brown, 2002; Bryant et al., 2012). It seems that water depth is only a proxy for multiple physiochemical variables, and we cannot completely exclude the influence of other environmental and biological factors on the vertical distribution of bacterial phyla.