Characterization of Bacterial Communities from the Surface and Adjacent Bottom Layers of Water in the Billings Reservoir

Here, we describe the bacterial diversity and physicochemical properties in freshwater samples from the surface and bottom layers of the Billings Reservoir, the largest open-air storage ecosystem in the São Paulo (Brazil) metropolitan area. Forty-four samples (22 from the surface and 22 from the bottom layers) were characterized based on 16S rRNA gene analysis using Illumina MiSeq. Taxonomical composition revealed an abundance of the Cyanobacteria phylum, followed by Proteobacteria, which were grouped into 1903 and 2689 different genera in the surface and the deep-water layers, respectively. Chroobacteria, Actinobacteria, Betaproteobacteria, and Alphaproteobacteria were the most dominant classes. The Shannon diversity index was in the range of 2.3–5.39 and 4.04–6.86 in the surface and bottom layers, respectively. Flavobacterium was the most predominant pathogenic genus. Temperature and phosphorus concentrations were among the most influential factors in shaping the microbial communities of both layers. Predictive functional analysis suggests that the reservoir is enriched in motility genes involved in flagellar assembly. The overall results provide new information on the diversity composition, ecological function, and health risks of the bacterial community detected in the Billings freshwater reservoir. The broad bacterial diversity indicates that the bacterioplankton communities in the reservoir were involved in multiple essential environmental processes.


Introduction
Access to safe, clean, and affordable water is one of the most basic humanitarian goals and is a major global challenge for the 21st century [1,2]. Both pollution from the outside and pollution from inside the reservoir worsen the water quality [3]. Exogenous contamination comes from sources such as industrial waste; sewage runoff; solid waste; meteoric water; and surface runoff from cities, farms, and forests [4][5][6]. Endogenous pollution comes from sediments, phytoplankton, aquatic plants, and aquaculture. Organic matter and biological and chemical oxygen demand (BOD and COD) in reservoirs can be either allochthonous or autochthonous. While the autochthonous organic matter is produced by phytoplankton and hydrophytes, allochthonous organic matter is introduced into aquatic systems by rainfall runoff [7]. BOD and COD levels indicate reservoir contamination, and rising nutrient levels cause reservoir BOD and COD levels to rise. Reservoirs are critical artificial aquatic bodies for many drinking water supply systems because they can maintain an equilibrium of storing or releasing water and play a central role in biogeochemical cycling, energy flows, and the recycling of nutrients. The functioning of these cycles is mainly maintained Life 2022, 12, 1280 2 of 17 by the inhabiting microorganisms. Despite their potential importance, the structure and functioning of these ecosystems' microbial communities have received less attention than other freshwater bodies, such as natural lakes and rivers [8].
Several studies have analyzed the structure and ecology of bacterial communities and have revealed the unexpected diversity of functional microbes in freshwater ecosystems [9][10][11][12][13][14][15][16]. For instance, a recent study by Shen et al. [17] has revealed the relationships between trophic status and planktonic microbial community structure and metabolic potential in freshwater lakes on the Yun-Gui Plateau, China. Wu K. et al. [18] showed that water depth was the most critical factor affecting the relative abundance of 11 dominant bacterial phyla, including Proteobacteria as the most dominant phylum in this lake, followed by Bacteroidetes, Chloroflexi, and Nitrospirae in the sediments of Lake Lugu in Southwest China. An earlier study by De Wever et al. [19] undertook a denaturing gradient gel electrophoresis analysis of PCR-amplified 16S RNA fragments approach to reveal the presence of pronounced vertical and latitudinal gradients in the bacterial community composition in Lake Tanganyika in East Africa. Actinobacteria and Gammaproteobacteria were found in the epilimnion, whereas Nitrospirae, green nonsulfur bacteria, Gammaproteobacteria, Actinobacteria, Deltaproteobacteria, and Firmicutes were found in the hypolimnion. The genetic diversity and activity of planktonic bacteria in freshwater systems are strongly shaped by the surrounding environmental conditions. Several studies have demonstrated that bacterial community structures depend on physical conditions and the availability of resources [20,21]. Moreover, it has been demonstrated that the extensive partitioning of water masses in such ecosystems may lead to progressively increased differences in both bacterial diversity [22] and water chemistry [23].
The Billings Reservoir in the Brazilian state of São Paulo is the largest freshwater storage aquatic body in the São Paulo Metropolitan Region, covering 127 km 2 . It has a total volume of 1228.7 × 106 m 3 and a maximum depth of 18 m [24,25]. The reservoir is used for a variety of purposes, including hydropower generation, water supply for 4.5 million people and industries, irrigation, fisheries, and flood control. The basin has a 20-km-long narrow central body that is connected and fed by eight branches, namely the Rio Grande, Rio Pequeno, Capivari, Pedra Branca, Taquacetuba, Bororé, Cocaia, and Alvarenga rivers [26][27][28]. Since the 1940s, the reservoir properties have changed significantly. At that time, part of the polluted water from the Tiete River (São Paulo city) was allowed to flow into the reservoir, believing that it would raise the water level to generate electric power [29]. This process, along with uncontrolled urban growth, significantly contributed to considerable anthropogenic eutrophication and an algal bloom in the reservoir [30].
There are very few studies on the bacterial community structures and compositions of the surface water in Brazil [31,32]. Thus, the goals of this study were to (i) explore the bacterial communities in surface and bottom layer water along the Billings Reservoir using 16S rRNA gene-based Illumina MiSeq sequencing, (ii) evaluate the presence of potential pathogens in these water samples, and (iii) explore the predicted functional biomarkers of the obtained bacterial communities in the basin to determine their role in the ecosystem.

Study Sites and Sample Collection
The study area covered the entire 127 km 2 of the Billings Reservoir, which is located west of the city of Sao Paulo at 23 • 47 S, 46 • 40 W, at an altitude of 746 m a.s.l. Surface and bottom layer water samples (250 mL each) were collected in July 2019 from 30 locations (approximately 17 km apart; Figure 1) in triplicate using a Van Dorn sampler as previously described [5]. Surface samples were collected from the side of a boat just below the water surface (10 cm) and the bottom water was collected at approximately 0.1 m above the sediment to avoid collecting soil material. Depending on the locations of the collected samples, the water depth of the study area from the reservoir ranged from 1.16 m to 13.6 m. A total of 180 samples were collected and transported to the laboratory on ice the same day they were collected. For the 16S rRNA gene analysis, duplicate samples from each Life 2022, 12, 1280 3 of 17 point were combined in 500 mL sterile bottles and thoroughly mixed. Then, a subsample of 40 mL was transferred into a sterile 50-mL tube (BD Falcon, Schaffhausen, Switzerland) and centrifuged for 10 min at 5000× g at 4 • C, and the supernatant was discarded, leaving approximately 600 µL of residual liquid. The pellets were frozen at −80 • C until DNA extraction. Water temperature (Temp), pH, dissolved oxygen (DO), and pH were measured on site from each sample using a handheld multiparameter water quality sonde (YSI Inc., Yellow Springs, OH, USA/Xylem Inc., Washington, DC, USA). Other variables, including turbidity, nitrate (NO 3 − ), sulfate (SO 4 2− ), phosphorus (P), and ammonia (NH 3 ), were determined from another subsample (250 mL) according to the Brazilian standard issued by the Environment National Council (CONAMA resolution 357/2005). Figure 1 provides the data of the sample sites.
(approximately 17 km apart; Figure 1) in triplicate using a Van Dorn sampler as previously described [5]. Surface samples were collected from the side of a boat just below the water surface (10 cm) and the bottom water was collected at approximately 0.1 m above the sediment to avoid collecting soil material. Depending on the locations of the collected samples, the water depth of the study area from the reservoir ranged from 1.16 m to 13.6 m. A total of 180 samples were collected and transported to the laboratory on ice the same day they were collected. For the 16S rRNA gene analysis, duplicate samples from each point were combined in 500 mL sterile bottles and thoroughly mixed. Then, a subsample of 40 mL was transferred into a sterile 50-mL tube (BD Falcon, Schaffhausen, Switzerland) and centrifuged for 10 min at 5000× g at 4 °C, and the supernatant was discarded, leaving approximately 600 µL of residual liquid. The pellets were frozen at −80 °C until DNA extraction. Water temperature (Temp), pH, dissolved oxygen (DO), and pH were measured on site from each sample using a handheld multiparameter water quality sonde (YSI Inc., Yellow Springs, OH, USA/Xylem Inc., Washington, DC, USA). Other variables, including turbidity, nitrate (NO3 − ), sulfate (SO4 2− ), phosphorus (P), and ammonia (NH3), were determined from another subsample (250 mL) according to the Brazilian standard issued by the Environment National Council (CONAMA resolution 357/2005). Figure 1 provides the data of the sample sites.

DNA Isolation and Library Preparation
Total genomic DNA was extracted from 200 µL resuspended pellets using the Pow-erSoil DNA kit (MO BIO Laboratories, Carlsbad, CA, USA), according to the manufacturer's instructions. To minimize potential bias during DNA extraction, each sample was extracted in duplicate and then pooled to quantify their DNA yield with a Qubit fluorometer (Invitrogen, Waltham, MA, USA). Sterile reagent grade type I water was used as a negative control for DNA extraction. The extracted DNA from each sample and the negative control was subjected to PCR amplification of the V3-V4 variable region of the 16S rRNA gene using the previously published primers Bakt_341F/Bakt_805R [33] and conditions previously described by our group [34,35]. A negative control containing all reagents except the template DNA was included with each set of reaction mixtures. Library preparation and massively parallel sequencing (MPS) were performed according to the manufacturer's protocol for paired-end sequencing using MiSeq Reagent Kit v3, as previously reported [5,34,35].

DNA Isolation and Library Preparation
Total genomic DNA was extracted from 200 µL resuspended pellets using the Power-Soil DNA kit (MO BIO Laboratories, Carlsbad, CA, USA), according to the manufacturer's instructions. To minimize potential bias during DNA extraction, each sample was extracted in duplicate and then pooled to quantify their DNA yield with a Qubit fluorometer (Invitrogen, Waltham, MA, USA). Sterile reagent grade type I water was used as a negative control for DNA extraction. The extracted DNA from each sample and the negative control was subjected to PCR amplification of the V3-V4 variable region of the 16S rRNA gene using the previously published primers Bakt_341F/Bakt_805R [33] and conditions previously described by our group [34,35]. A negative control containing all reagents except the template DNA was included with each set of reaction mixtures. Library preparation and massively parallel sequencing (MPS) were performed according to the manufacturer's protocol for paired-end sequencing using MiSeq Reagent Kit v3, as previously reported [5,34,35]. ATT TGT GTA GAA GGT GC-3 and mcyE-R4; 5 -AAT TCT AAA GCC CAA AGA CG-3) [37] were used. Amplification of these fragments was conducted using a 50-100 ng DNA template, 2 mM MgCl2, 0.1 mM dNTPs, 0.5 µM of each primer, and 2.5 U high-fidelity Taq platinum DNA polymerase (Invitrogen, Carlsbad, CA, USA) in a MgSO 4 reaction buffer.

Detection of Toxin-Producing Cyanobacterial Genes
After an initial denaturation of 5 min at 94 • C, 35 cycles of 30 s at 94 • C, 30 s at 55 • C, 60 s at 72 • C, and a final extension at 72 • C for 5 min were performed. Each PCR included a known cyanobacterial DNA positive control and an interspersed no DNA template negative control. The amplified product was electrophoresed through 1% (wt/vol) agarose gels containing 0.5× Tris Borate EDTA, followed by ethidium bromide staining.

Bioinformatics and Statistical Analysis
Base-calling and data quality were initially assessed on an MiSeq instrument using RTA v1.18.54 and MiSeq Reporter v2.6.2.3 software (Illumina Inc., San Diego, CA, USA). The sequences were analyzed by a pipeline of the 16S Microbiome Taxonomic Profiling (MTP) of the EzBioCloud (https://www.ezbiocloud.net/ (accessed on 19 March 2020)) application update 2019.04.09. Briefly, the processing of raw reads started with quality checking and filtering of low quality (<Q25) reads by Trimmomatic version 0.32 [38]. After QC pass, paired-end sequence data were merged together using the fastq_mergepairs command of VSEARCH version 2.13.4 [39] with default parameters. Primers were then trimmed with the alignment algorithm of Myers and Miller [40] at a similarity cut off of 0.8. Nonspecific amplicons that do not encode 16S rRNA were detected by Nhmmer [41] in the HMMER software package version 3.2.1 with hmm profiles. Unique reads were extracted, and redundant reads were clustered with the unique reads by the derep_fulllength command of VSEARCH [39]. The EzBioCloud 16S rRNA database [42] was used for taxonomic assignment using the usearch_global command of VSEARCH [39], followed by more precise pairwise alignment [40]. Chimeric reads were filtered on reads with <97% similarity by reference-based chimeric detection using the UCHIME algorithm [43] and the nonchimeric 16S rRNA database from EzBioCloud. Following chimeric filtering, reads that could not be identified to the species level (with <97% similarity) in the EzBioCloud database were compiled, and the cluster_fast command [39] was used to perform de novo clustering to generate additional OTUs. Finally, OTUs with single reads (singletons) were omitted from further analysis. All analyses of molecular data were performed on both nonrarefied and rarefied data to 15,766 sequences (the size of the smallest library). Our initial analysis revealed that the most predominant bacterioplankton taxa and their abundances in the rarefied data were very similar to those in the nonrarefied data. To avoid discarding data, we decided to run our analysis on nonrarefied data. The ACE [44], Chao1 [45], and jackknife [46] α-diversity indices were used to calculate the bacterial richness, and the Shannon [47], Simpson function [47], and NPShannon [48] indices were used to estimate the bacterial evenness in each group using the Wilcoxon rank-sum test. β-diversity was calculated using the Bray-Curtis dissimilarity distance [49] to reveal differences between surface and bottom layer bacteriomes based on their OTU data. The predicted profiles were categorized into clusters of Kyoto Encyclopedia of Genes and Genomes (KEGG) orthology and pathways. To determine whether the reservoir bacterial compositions are affected by depth status, we performed a principal coordinate analysis (PCoA) at the genus level using the Bray-Curtis dissimilarity distance. A pairwise permutational multivariate analysis of variance (PERMANOVA) test was computed to determine whether there were differences in the bacterial communities between the two layers. The enrichment in the assigned taxonomic and functional profiles of the two groups was defined by linear discriminant analysis (LDA) of the effect size (LEfSe) algorithm. All of the above analyses were carried out in ChunLab's bioinformatics cloud platform, EzBioCloud 16S-based MTP. Statistically significant differences in physicochemical characteristics between both layers were determined by using Student's t tests and visualized in GraphPad Prism software V8. The correlations between bacterial community diversity (richness, evenness, and OTUs) and water properties were assessed using bootstrap-based confidence limits (100 iterations) and  [50]. Uncorrected and corrected p values of ≤0.05 were used for all statistical analyses.
For the detection of bacterial pathogens, we considered any bacteria potentially pathogenic if at least one species with a minimum abundance of 10 strains of any genus was categorized as biosafety level 2 or 3 by the American Biological Safety Association (https://my.absa.org (accessed on 7 May 2020)). The analyses were performed based on uncultivated microorganisms.

Physicochemical Characteristics of Water Samples
The physicochemical properties of water samples collected from the surface and bottom layers of the Billings Reservoir at 30 different locations are demonstrated in Tables 1 and 2. The samples were collected in July 2019 when the temperature was above 20 • C and the season was mildly dry. The temperature from both layers ranged from 18.8 • C to 22.1 • C, and DO ranged from 3.5 to 9.5 mg/L. The higher DO levels at some sampling locations are most likely due to wind and aeration. Both temperature and DO parameters were significantly higher (p < 0.05) in surface water than in the bottom layers (average temperature: 21.1 • C vs. 20.7 • C, averaged DO: 9.5 mg/L vs. 8.4 mg/L), whereas phosphorus concentrations showed the opposite pattern. The concentrations of ammonia and nitrate along the reservoir showed no obvious difference between the surface and bottom layers. PC1 explained a large proportion of the total variance (29.6%), and the proportion increased to 46.3% when PC2 was added ( Figure 2). The variables contributing most to the PC1 were CHAO (r = 0.50) and NO 3 − (r = 0.41), while those contributing most to the PC2 were DO (r = −0.6) and PH (r = −0.52). The bacterial diversity and OTU relative abundance in the Billings Reservoir were positively correlated with the NO 3− concentration.

Bacterial Community Structure
Of the 30 locations, 22 paired water samples from the surface and bottom layers were successfully amplified, sequenced, and submitted for 16S rRNA analysis. Up to 100,000 MPS reads from each sample were uploaded, quality controlled, and profiled by the EzBi-oCloud tool. The trimming-based quality control removed 30 Tables S1 and S2). Rarefaction curves reached stable values, indicating that the diversity of bacterial populations in both groups was sufficiently covered by the generated sequences (data not shown). The number of OTUs per sample ranged from 1063 to 5919 and from 305 to 2150 in the surface and bottom layer samples, respectively. The α diversity of bacteria in the surface and bottom layers of the reservoir was compared by both richness (i.e., the observed OTUs, ACE,   Tables S1 and S2). Rarefaction curves reached stable values, indicating that the diversity of bacterial populations in both groups was sufficiently covered by the generated sequences (data not shown). The number of OTUs per sample ranged from 1063 to 5919 and from 305 to 2150 in the surface and bottom layer samples, respectively. The α diversity of bacteria in the surface and bottom layers of the reservoir was compared by both richness (i.e., the observed OTUs, ACE, Chao1, and jackknife index) and evenness indices (i.e., the NPShannon, Shannon, Simpson, and phylogenetic diversity index). The results showed that, with the exception of the phylogenetic diversity index, no α indices differed significantly between the bacterial communities from both layers. The phylogenetic diversity index was significantly higher in the bottom layers than in the surface layers (Wilcoxon rank-sum test; Chao1, and jackknife index) and evenness indices (i.e., the NPShannon, Shannon, Simpson, and phylogenetic diversity index). The results showed that, with the exception of the phylogenetic diversity index, no α indices differed significantly between the bacterial communities from both layers. The phylogenetic diversity index was significantly higher in the bottom layers than in the surface layers (Wilcoxon rank-sum test; p < 0.005) ( Figure  3). The principal coordinate analysis of the β-diversity d results revealed no significant (Wilcoxon rank-sum test; p = 0.131) similarities between the bacterial communities from the reservoir's surface and bottom layers.

Identification of the Billings Bacteriome in Surface and Bottom Layer Samples
Seven phyla (Cyanobacteria, Proteobacteria, Actinobacteria, Bacteroidetes, Verrucomicrobia, Planctomycetes, and Chlorobi) were found to be abundant in samples from both layers.  (Supplementary Tables S3 and S4). At the order level, 21.9% and 16.9% of OTUs were assigned to Oscillatoriales in the surface and bottom layers, respectively. The order Planktophila accounted for >12% of all detected OTUs from Billings cyanobacterial sources, followed by the order Burkholderiales, which accounted for >24%. At the genus level, 20.2% and 14.5% of OTUs were assigned to the Planktothrix genus within the Microcoleaceae family in the surface and bottom layers, respectively, and >8% OTUs were assigned to the Nanopelagicus genus within the family Nanopelagicaceae in both layers ( Figure  5). Among the Cyanobacteria, Chroobacteria were the most abundant group in the reservoir, with a median relative abundance of 18.96% and 18.08% of all OTUs in the surface and bottom layer samples, respectively. Differentially abundant taxa between the two water samples were identified using the LEfSe algorithm with a minimum LDA score of 2.0 [51]. This analysis revealed 126 taxa, which included 11 classes, 24 families, 29 genera, 24 orders, 6 phyla, and 32 species. Of these, two families, four genera, one order, and nine

Identification of the Billings Bacteriome in Surface and Bottom Layer Samples
Seven phyla (Cyanobacteria, Proteobacteria, Actinobacteria, Bacteroidetes, Verrucomicrobia, Planctomycetes, and Chlorobi) were found to be abundant in samples from both layers. The relative abundance of taxa fluctuated among surface and bottom layer samples at the phylum level, as shown in Supplementary Figures S1 and S2. Cyanobacteria predominated the bacteriome of the water in the Billings Reservoir, but the individual composition datasets showed a significant amount of variation. The proportion of cyanobacteria ranged from 1.4% (B04S) to 71.1% (B09S) in the surface layer, and from 0.4% (B04F) to 44.6% (B09S) in the bottom layer. The proportion of the phylum Proteobacteria ranged from 10.2% (B09S) to 47% (B04S) in the surface layers and from 18.9% (B09S) to 56% (B04S) in the bottom layers. Proteobacteria and Bacteroidetes were less abundant in the surface layers compared to the bottom layer of the samples (Figure 4). The most common classes in all samples were Chroobacteria, Actinobacteria, Betaproteobacteria, Alphaproteobacteria, Sphingobacteria, and Acidimicrobiia (Supplementary Tables S3 and S4). At the order level, 21.9% and 16.9% of OTUs were assigned to Oscillatoriales in the surface and bottom layers, respectively. The order Planktophila accounted for >12% of all detected OTUs from Billings cyanobacterial sources, followed by the order Burkholderiales, which accounted for >24%. At the genus level, 20.2% and 14.5% of OTUs were assigned to the Planktothrix genus within the Microcoleaceae family in the surface and bottom layers, respectively, and >8% OTUs were assigned to the Nanopelagicus genus within the family Nanopelagicaceae in both layers ( Figure 5). Among the Cyanobacteria, Chroobacteria were the most abundant group in the reservoir, with a median relative abundance of 18.96% and 18.08% of all OTUs in the surface and bottom layer samples, respectively. Differentially abundant taxa between the two water samples were identified using the LEfSe algorithm with a minimum LDA score of 2.0 [51]. This analysis revealed 126 taxa, which included 11 classes, 24 families, 29 genera, 24 orders, 6 phyla, and 32 species. Of these, two families, four genera, one order, and nine species were significantly abundant and discriminative between the groups (FDR-adjusted p value < 0.02. Supplementary Table S5). The relative abundance of bacteria from the genera Stenotrophomonas, Achromobacter, Comamonas, Pseudomonas, Acinetobacter, and Schlesneria was significantly higher in the bottom layers than in the surface samples (FDR-adjusted p-value < 0.05). species were significantly abundant and discriminative between the groups (FDR-adjusted p value < 0.02. Supplementary Table S5). The relative abundance of bacteria from the genera Stenotrophomonas, Achromobacter, Comamonas, Pseudomonas, Acinetobacter, and Schlesneria was significantly higher in the bottom layers than in the surface samples (FDRadjusted p-value < 0.05).

Prediction of Functional Biomarkers
An LEfSe analysis was conducted to predict the most pertinent functional pathways that shape the Billings bacteriome between both layers. The software package PICRUSt2, implemented in the EzBiocloud online tool, was used to infer the content of bacterial genes from the 16S rRNA data and aggregate the relative abundance of functional genes into metabolic pathways. From all OTUs detected in matched samples, 83 KEGG orthology (KO) terms were predicted. Of these, four differentially abundant (FDR-adjusted p < 0.05) KO terms between surface and bottom layer water were identified (Supplementary Table  S6), and all had LDA effect sizes <3. Most of the differentially abundant predicted KO terms, including protein metabolism, signaling and cellular processes, and cell motility were most abundant in bottom layer samples. The pathway of the flagellar assembly was also enriched in the bottom water layer group, based on the cutoff of uncorrected p < 0.02 instead of FDR corrected <0.89, implying that the growth environment for the bacterial communities in the bottom layer water was much higher than that of the surface layer water (LDA = 2.24). PICRUSt module analysis demonstrated that the microbiota in the bottom layers exhibits increased biosynthesis of tetrahydrofolate (M00841) and C21- species were significantly abundant and discriminative between the groups (FDR-adjusted p value < 0.02. Supplementary Table S5). The relative abundance of bacteria from the genera Stenotrophomonas, Achromobacter, Comamonas, Pseudomonas, Acinetobacter, and Schlesneria was significantly higher in the bottom layers than in the surface samples (FDRadjusted p-value < 0.05).

Prediction of Functional Biomarkers
An LEfSe analysis was conducted to predict the most pertinent functional pathways that shape the Billings bacteriome between both layers. The software package PICRUSt2, implemented in the EzBiocloud online tool, was used to infer the content of bacterial genes from the 16S rRNA data and aggregate the relative abundance of functional genes into metabolic pathways. From all OTUs detected in matched samples, 83 KEGG orthology (KO) terms were predicted. Of these, four differentially abundant (FDR-adjusted p < 0.05) KO terms between surface and bottom layer water were identified (Supplementary Table  S6), and all had LDA effect sizes <3. Most of the differentially abundant predicted KO terms, including protein metabolism, signaling and cellular processes, and cell motility were most abundant in bottom layer samples. The pathway of the flagellar assembly was also enriched in the bottom water layer group, based on the cutoff of uncorrected p < 0.02 instead of FDR corrected <0.89, implying that the growth environment for the bacterial communities in the bottom layer water was much higher than that of the surface layer water (LDA = 2.24). PICRUSt module analysis demonstrated that the microbiota in the bottom layers exhibits increased biosynthesis of tetrahydrofolate (M00841) and C21-

Prediction of Functional Biomarkers
An LEfSe analysis was conducted to predict the most pertinent functional pathways that shape the Billings bacteriome between both layers. The software package PICRUSt2, implemented in the EzBiocloud online tool, was used to infer the content of bacterial genes from the 16S rRNA data and aggregate the relative abundance of functional genes into metabolic pathways. From all OTUs detected in matched samples, 83 KEGG orthology (KO) terms were predicted. Of these, four differentially abundant (FDR-adjusted p < 0.05) KO terms between surface and bottom layer water were identified (Supplementary Table S6), and all had LDA effect sizes <3. Most of the differentially abundant predicted KO terms, including protein metabolism, signaling and cellular processes, and cell motility were most abundant in bottom layer samples. The pathway of the flagellar assembly was also enriched in the bottom water layer group, based on the cutoff of uncorrected p < 0.02 instead of FDR corrected <0.89, implying that the growth environment for the bacterial communities in the bottom layer water was much higher than that of the surface layer water (LDA = 2.24). PICRUSt module analysis demonstrated that the microbiota in the bottom layers exhibits increased biosynthesis of tetrahydrofolate (M00841) and C21-steroid hormone (M00109), while the surface microbiota showed an increased use of the Mce transport system (M00670).

Search for Predefined Bacterial Groups, Pathogens, and Cyanotoxin Genes
Screening for pre-defined bacterial groups in the surface and bottom layers of the Billings Reservoir revealed important taxa associated with the human gut, including the phylum Proteobacteria (surface; median relative abundance value, 24.7%, bottom layer; median relative abundance value, 27.7%). The search for bacterial pathogens in both layers identified several bacterial genera known to be human, animal, and/or plant pathogens, including Flavobacterium, Legionella, Staphylococcus, Pseudomonas, Aeromonas, Acidovorax, and Acinetobacter. Flavobacterium (surface; median relative abundance value, 0.65%, bottom layer; median relative abundance value, 0.94%) and Legionella (surface; median relative abundance value, 0.21%, bottom layers; median relative abundance value, 0.19%) were the predominant genera.
The DNA of the three selected genes involved in the biosynthesis of mcy (mcyA, mcyE, and mcyD) was successfully amplified in all analyzed samples, confirming the genetic potential of the strains in the reservoir to produce microcystin (data not shown).

Distribution and Diversities of Bacterial Communities along the Billings Reservoir
The bacterial diversity in surface water from the Billings Reservoir has recently been studied, with an emphasis on bacterial structure; nonetheless, detailed local bacterial heterogeneity remains unknown. In this study, we investigated the distribution of bacterial communities in the surface and bottom layers of the Billings Reservoir, as well as the environmental conditions that influence their composition. Our results demonstrated that the overall α bacterial diversity represented by the Shannon index was relatively higher in the Billings Reservoir than in other studied freshwater systems: Itupararanga Reservoir, State of São Paulo, Brazil [52], in Pavin Lake [53], and Mayinghai, Pipahai, and Gonghai Lakes in China [54]. Perhaps the greatest value of biodiversity is attributed to intensive anthropogenic disturbance in the reservoir, which may have resulted in increased nutrient discharges, resulting in greater nitrogen and phosphorus concentrations, as indicated in this study. Dias et al. [55] demonstrated that the influence of surrounding human activities, as measured by the human footprint index, had a significant impact on the functional diversity and trait composition of fish assemblages in the Billings Reservoir. It is possible that other variables not included in our study could have contributed to the variability of the bacterial community in the Billings Reservoir.
The observed bacterial community structure and diversity indices in this study were comparable to those in other reservoirs [56]. The taxonomical composition revealed that the Cyanobacteria phylum was the most abundant in the community, followed by Proteobacteria and Actinobacteria. These results are consistent with other Brazilian studies that describe bacterial communities in the Amazon basin [57][58][59] and Tocantins River [60], but are somewhat different from those previously described for reservoirs in China. For instance, Qu J et al. [56] used the 16S rDNA Illumina approach to show that Firmicutes, Proteobacteria, Cyanobacteria, and Bacteroidetes were the dominant bacterial phyla in the Miyun Reservoir, which is considered the largest man-made reservoir in North China. These differences in dominating phyla may be due to a variety of environmental factors, such as air, soil, and water pollution, rainfall-induced nutrient fluctuations, and changes in local conditions. The difference may also be attributed to the use of distinct DNA extraction methods and/or primer selection [61] or different 16S rDNA gene bioinformatic pipelines. Although the bacterial richness indices varied little between the surface and bottom layer samples in this study, some specific bacterial groups showed a clear difference. For instance, members of the β-Proteobacteria and γ-proteobacteria, Acidobacteria, Bacteroidetes, and Sphingobacteriales were significantly more abundant in the bottom layer than in the surface water. Because the water flow within certain layers could contribute significantly to the displacement of bacterial populations, it is reasonable to assume that the latter could constitute a contribution to the observed bacterial community differences between the layers. It is also conceivable that the phosphorus availability in the bottom layers compared to the surface water could positively influence the growth rates of these bacterial groups [62].
In addition to the seven dominant phyla in the Billings Reservoir, there were approximately 3.5% and 4.4% of unclassified taxa in the surface and bottom layers water, respectively, indicating that as-yet-unidentified bacterial populations with unknown metabolic

Dominant Taxonomic Groups in the Bacterioplankton Community of the Reservoir
Members of Cyanobacteria were detected as the most dominant genus, and all Cyanobacteria blooms tested here were toxic. These results lend further support to previous studies that found cyanobacterial communities in Brazilian semiarid reservoirs [32,[63][64][65][66][67]. The presence of these bacteria in high relative abundance in the reservoir can be linked to uploading nutrients such as ammonia and an increase in water temperature (20.7 • C) at the time of sampling. The OTUs of these bacteria have been shown to outcompete other planktonic microbes for nutrients in eutrophic systems [68]. A previous study on the Neuse River, North Carolina, conducted by Paerl, recommended that a reduction of 30-40% of NO 3 had the optimum power of minimizing M. aeruginosa as a dominant phytoplankter [69]. Concerning temperature, there is a consensus among researchers that water temperatures below 20 • C are generally considered unfavorable for the development of common water bloom-forming genera such as Dolichospermum and Microcystis. In contrast to this, our findings showed widespread Microcystis and Planktothrix in the reservoir. The presence of more Cyanobacteria in freshwater systems has a variety of impacts, including oxygen depletion, fish mortality, and toxicity, which offer a variety of health risks [70]. Because the Cyanobacteria detected in this study are capable of producing toxins, their presence in the reservoir poses a huge threat to animal and human health.
In this study, Proteobacteria constituted a sizable fraction of the bacterial community in the reservoir, which is consistent with previous reports [71,72]. Despite the dominance of Proteobacteria in both layers, a significant difference in the composition of these members was observed, as detailed in the taxonomic analysis. Betaproteobacteria were the most frequently detected Proteobacteria group in this work, which agrees with other studies [73,74]. Microorganisms belonging to the Betaproteobacteria class have been associated with anthropogenic activities [75]. The alpha and Gammaproteobacteria found in both layers probably indicate an increase in organic and inorganic inputs, as well as phytoplankton production [62,73].
The Actinobacteria phylum was the reservoir's third most abundant OTU phylum. Members of this phylum are among the most abundant groups in freshwater habitats [73,76]. The relative abundance of these microbes is inversely correlated to that of Cyanobacteria, which cause prolonged and irretrievable ecological disturbances to freshwater ecosystems and serve as sentinels of impending ecological damage [77,78]. The OTUs of Actinobacteria were less abundant than those of Proteobacteria, which comprised 23.7% and 28.5% of the total bacterial relative abundance in the surface and bottom layer water of the Billings Reservoir, respectively. Actinobacteria, which are well-recognized soil bacteria [79], are frequently detected in oligotrophic freshwater habitats [80] and are often associated with oligotrophic ecosystems [81]. They have long been known to produce pigments that protect them from UV radiation, which may easily penetrate deep into a freshwater habitat [73].

Bacterioplankton Community Diversity between the Surface and Bottom Layers
The alpha diversity analysis conducted between the surface and bottom layers revealed high phylogenetic diversity indices of the bacteriome inhabiting the bottom layers. This result may indicate that bacterial communities in the bottom layers have experienced high diversification rates or the successful immigration of multiple lineages. One factor that may promote higher bacterial diversification in the bottom layers is that this habitat is possibly less extreme than that of the surface layers and permits easier radiation [82]. An examination of Billings water did not reveal significant differences in bacterial community β-diversity between the two layers, suggesting that a core bacteriome exists between both layers of the reservoir.

Potential Functions and Pathways of Bacterial Communities in the Reservoir
Application of the LEfSe method demonstrated the presence of Stenotrophomonas species as the most significant specific biomarkers in the bottom layers. These species play an important ecological role in the nitrogen and sulfur cycles, and several Stenotrophomonas species can engage in beneficial interactions with plants, promoting growth and protecting plants from attack [83]. Because the reservoir contains a large number of Cyanobacteria, the death of many of these bacteria could result in the release of a high contents of sulfur-containing amino acids from their cells [84], resulting in sulfur-rich water and thus explaining the detection of Stenotrophomonas species. However, these conclusions must be interpreted with caution given the serious limitations of deducing the biofunctions of the reservoir bacteria from 16S rRNA gene sequencing and the potential risk of bias.
In addition to taxonomic compositions, we identified a pathway containing genes encoding for flagellar assembly that were abundant in the Billings 16S metagenome. The occurrence of flagellar assembly genes within the bottom layers indicates that the reservoir contains bacterial communities that utilize flagella for locomotion. It is possible that the adhesion process of these bacteria was influenced by environmental factors such as pH and higher levels of metals introduced into the reservoir because of the rapid expansion of industry and intense domestic activities.

Occurrence of Pathogenic Bacteria
Many developing countries are suffering widespread contamination of their freshwater supplies by bacterial pathogens, which have been linked to epidemics of numerous waterborne diseases [85]. In the current study, we found that pathogenic bacteria were ubiquitous across all the sampled waters in the Billings Reservoir. Among these pathogens, the Flavobacterium and Legionella genera were the most prevalent in both layers. Although most Flavobacteria are harmless, some are opportunistic or true pathogens that cause diseases in animals, plants, and humans [86][87][88]. The majority of Legionella species are considered pathogenic and have been reported as one of the leading bacterial etiological pathogens of waterborne outbreaks in the United States between 2007 and 2009 [89].

Limitations and Conclusions
The main limitation of this study is that we restricted the analysis to samples at a single point in time because water samples from the bottom layers were collected only once. While the spatiotemporal data from surface water are still incomplete, the analysis of the data that is now available is well underway. In addition, we used bacterial DNA genomics for this investigation, which would have revealed the presence of bacterial populations regardless of whether they were dead or alive, culturable cells, or nonculturable cells. Another limitation, inherent in 16S rRNA diversity surveillance by the Illumina MiSeq platform, is the fact that the taxonomic and ecological resolution of the 16S rRNA gene V3-V4 region is insufficient for accurately identifying pathogenic bacteria and strains with flagella genes to species levels, as revealed in functional analysis. These results should be amenable to robust metagenomic data, RNA sequencing, and cultivable fractions to expand our understanding of the bacterial strains associated with the reservoir, their genes, and metabolic pathways. That being said, our findings do reveal the extent of bacterial diversity in the reservoir, the health hazards for household drinking water, and the effects of the direct use of water by residents of the reservoir communities.
To summarize, the toxin-producing Cyanobacteria phylum was found to be the most abundant in the community, followed by Proteobacteria and Actinobacteria. Ammonia and high water temperature are two potential contributors to the abundance of Cyanobacteria in the reservoir. According to functional prediction analyses, genes associated with cell motility, such as flagellar assembly, were abundant in the Billings bacteriome.
In conclusion, this study provides important information about the numerous bacteria inhabiting the Billings Reservoir and the combination of environmental factors that shape their structure. These results may help pave the way for future studies devoted to control-ling and improving the water quality in the Billings Reservoir, which is facing rapid urban development and urbanization.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/life12081280/s1, Figure S1: The relative abundance of dominant taxa at the phylum level in surface water samples. Figure S2: The relative abundance of dominant taxa at the phylum level in bottom-layer samples. Table S1: Valid reads sequenced for each matched sample from the surface layer and alpha diversity indices at the 97% OTU level of 16S rRNA gene fragments. Table S2: Valid reads sequenced for each matched sample from the bottom layer and alpha diversity indices at the 97% OTU level of 16S rRNA gene fragments. Table S3: Count of the individual surface layer bacterial microbiota of Billings reservoir on a phylum, class, order, family, and genus levels. Table S4: Count of the individual bottom layer bacterial microbiota of Billings reservoir on a phylum, class, order, family, and genus levels. Table S5: Differentially abundant taxa between surface and bottom layer water samples identified using LEfSe analysis (minimum LDA score: 2.0).