Vertical Distribution of Bacterial Community in Water Columns of Reservoirs With Different Trophic Conditions During Thermal Stratification

Water eutrophication is a global ecological issue, and thermal stratification of water bodies can enable eutrophication. We examined bacterial communities in the stratified water columns and sediments in two different trophic reservoirs along the Wujiang River using quantitative real-time PCR and high-throughput sequencing. Bacterial 16S rRNA gene copies varied from 3.70 × 107 to 5.27 × 108 copies/L in the water column of Hongjiadu (HJD) Reservoir (60 m water depth) with slightly stratified variation; while in Wujiangdu (WJD) Reservoir (70 m water depth), bacterial abundance decreased markedly from the surface to the bottom(1.74 × 109 to 2.38 × 107 copies/L). The vertical distribution patterns of bacteria in both reservoirs resembled those of water Chlorophyll a (Chla) concentrations. The abundance was negatively correlated with water depth (D), total nitrogen (TN), nitrate (NO3–-N), and positively correlated with water temperature (T) and dissolved oxygen (DO) level. In contrast, the alpha diversity of bacteria showed the opposite trend in the vertical water column. Proteobacteria, Actinobacteria, and Bacteroidetes were the predominant phyla in the water column of both reservoirs. Compared to WJD Reservoir, HJD Reservoir displayed marked vertical spatial difference in bacterial community structure during thermal stratification. In particular, Pseudomonas was frequently detected at the bottom of the HJD Reservoir. These results were consistent with predictive metagenomic profiling that revealed different vertical functional variation patterns of the bacterial communities in the two reservoirs. The bacterial community structure of HJD Reservoir was associated with water D, ammonium (NH4+-N), nitrite (NO2–-N), and total phosphorus (TP). The community structure of WJD Reservoir was related to water T, Chla, NO3–-N, and TN. The findings highlighted the important roles played by thermal stratification and nutrients in shaping the water bacterial community structure. Additionally, the absolute abundance of water nitrifiers (AOB gene copies) and denitrifiers (narG, nirS, norB, and nosZ gene copies) displayed significant vertical differences in the water columns of both reservoirs. Gene copies involved in denitrification were significantly higher than those involved in nitrification. Water phosphorus and nitrogen contents were important variables influencing the absolute abundance of ammonia oxidizers and denitrifying bacteria, respectively. Our study revealed that the emergence of thermal stratification was responsible for the vertical stratification of bacteria in water and affected the bacterial community structure together with nutrients.

Water eutrophication is a global ecological issue, and thermal stratification of water bodies can enable eutrophication. We examined bacterial communities in the stratified water columns and sediments in two different trophic reservoirs along the Wujiang River using quantitative real-time PCR and high-throughput sequencing. Bacterial 16S rRNA gene copies varied from 3.70 × 10 7 to 5.27 × 10 8 copies/L in the water column of Hongjiadu (HJD) Reservoir (60 m water depth) with slightly stratified variation; while in Wujiangdu (WJD) Reservoir (70 m water depth), bacterial abundance decreased markedly from the surface to the bottom(1.74 × 10 9 to 2.38 × 10 7 copies/L). The vertical distribution patterns of bacteria in both reservoirs resembled those of water Chlorophyll a (Chla) concentrations. The abundance was negatively correlated with water depth (D), total nitrogen (TN), nitrate (NO 3 − -N), and positively correlated with water temperature (T) and dissolved oxygen (DO) level. In contrast, the alpha diversity of bacteria showed the opposite trend in the vertical water column. Proteobacteria, Actinobacteria, and Bacteroidetes were the predominant phyla in the water column of both reservoirs. Compared to WJD Reservoir, HJD Reservoir displayed marked vertical spatial difference in bacterial community structure during thermal stratification. In particular, Pseudomonas was frequently detected at the bottom of the HJD Reservoir. These results were consistent with predictive metagenomic profiling that revealed different vertical functional variation patterns of the bacterial communities in the two reservoirs. The bacterial community structure of HJD Reservoir was associated with water D, ammonium (NH 4 + -N), nitrite (NO 2 − -N), and total phosphorus (TP). The community structure of WJD Reservoir was related to water T, Chla, NO 3 − -N, and TN. The findings highlighted the important roles played by thermal stratification and nutrients in shaping the water bacterial community structure. Additionally, the absolute abundance of water nitrifiers (AOB gene copies) and denitrifiers (narG, nirS, norB, and nosZ gene copies) displayed significant vertical differences in the water columns of both reservoirs. Gene copies involved in denitrification were significantly higher than those involved in nitrification. Water phosphorus and nitrogen contents were important

INTRODUCTION
Water eutrophication has become one of the major environmental issues affecting lakes and reservoirs on a global scale (Paerl et al., 2016;Wang et al., 2019). Compared to natural rivers, artificial reservoirs generally exhibit hydrodynamic variations and longer exchange periods (Jardim et al., 2020). Additionally, numerous nutrients, such as nitrogen (N) and phosphorus (P) are intercepted, which alter the nutrient cycle, drive the rapid propagation of algae and other plankton, and decrease the dissolved oxygen (DO) content in the water body (Van Cappellen and Maavara, 2016;Shi et al., 2017;Petty et al., 2020), thereby worsening the eutrophication of the water body. According to 2018 data from the United Nations Water Resources Organization, over 75% of closed water bodies (such as lakes and reservoirs) are eutrophic to a certain extent, with deterioration in water quality.
Thermal stratification is a common phenomenon after reservoir completion. Because of the exchange of material and energy by the atmosphere and solar radiation, the temperature difference between the upper and lower layers of the water body leads to a difference in water density, which results in the formation of an epilimnion -thermocline -hypolimnion vertical profile from the top to the bottom of the water column (Tuo et al., 2014;Kirillin and Shatwell, 2016). The thermocline hinders the mixing of upper and lower waters and slows the vertical transfer of DO and particulates in the epilimnion layer , which causes anoxia in the bottom water and the release of nutrients from sediments, and aggravates the eutrophication (Conley et al., 2009;Rigosi and Rueda, 2012;Müller et al., 2016). With the appearance of thermal stratification, the vertical changes in pH, DO, turbidity, and Chla concentration are highly consistent and display marked stratification. Water stratification may be one of the major factors determining the vertical distribution of phytoplankton, which leads to the death of aquatic organisms, such as fish, in affected lakes. However, little is known about the effect of water thermal stratification on the abundance of bacteria in reservoirs (Yu et al., 2014;Wang et al., 2015;Yang et al., 2015).
In aquatic ecosystems, microorganisms are highly sensitive to changes in the water environment and are affected by a variety of environmental factors, including organic matter, nitrogen, phosphorus, pH, DO, and temperature (Zhao et al., 2012;Daggett et al., 2015;Ruiz-González et al., 2015;Sinistro et al., 2015;Li et al., 2017). Moreover, microbial activity drives nutrient cycling (Salcher et al., 2010) and influences water quality (Houshmand et al., 2020) to maintain the structure and function of aquatic ecosystems. Nitrogen serves as a limiting factor for water eutrophication (Canfield et al., 2010). Denitrification plays a critical role in nitrogen removal (Peng et al., 2016), and denitrifying bacteria contribute to over 87% of the nitrogen loss in freshwater (Schubert et al., 2006;Annika and Christopher, 2010;Chen et al., 2012). Denitrifiers are facultative anaerobes that are affected by DO (Oshiki et al., 2016), organic matter (Cornwell et al., 2014), and N availability (Smith et al., 2015). Examining the expression of the genes encoding for denitrifiers (such as NAR/NAP, NIR, NOR, and NOS) in the bacterial communities can help clarify the denitrification potential in freshwater ecosystems (Braker et al., 2015;Smith et al., 2015;Tatti et al., 2015;Gao et al., 2016;Mao et al., 2017).
In this study, we focused on the bacterial communities in two reservoirs of the Wujiang River in Southwestern China. Water columns and surface sediments from the lacustrine zone were chosen for comparison. The two major objectives were to (1) evaluate the effects of water thermal stratification on bacterial communities with various levels of eutrophication, and to (2) investigate the vertical dynamic distribution of nitrification and denitrification in the water columns.

Site Description and Sampling
The sampling areas of HJD and WJD reservoirs were located in the Wujiang River. The Wujiang River is the largest tributary south of the upper reaches of the Yangtze River. HJD Reservoir was built in 2004, with a total reservoir capacity of 49.3 × 10 8 m 3 , the highest water level of 1140 m, the normal water level area of 80.5 km 2 , and the retention time of the water body is 360 days, which suggests it is a regulating reservoir in the year of oligotrophic nutrition (Zhao et al., 2017); WJD Reservoir was built in 1979, with a drainage area of 27,790 km 2 , a total reservoir capacity of 21.4 × 10 8 m 3 , and the normal water level of 47.8 km 2 . The reservoir is a seasonal regulating reservoir with eutrophic nutrition (Zhao et al., 2017) and appreciably affected by human development and use. The samples were collected from the lacustrine zones of the two reservoirs of the Wujiang River in April 2017.
Water samples from the vertical profiles in both reservoirs (HJD and WJD) were collected at a range of depths (within 0,5,15,30,and 50 m for HJD;within 0,5,15,30,45, and 60 m for WJD) using a Niskin water sampler (Model 1010, General Oceanics, United States) (Figure 1). About 500 mL sample volumes were filtered through a 0.22 µm sterilized filter membrane (MF-Millipore, United States) and stored in 50 mL sterile centrifuge tubes at −80 • C. The remaining water samples were stored at 4 • C for physicochemical analyses. Sediment samples (top 0-20 cm) were collected using a Van Veen grab sampler (35 cm × 25 cm × 30 cm), and three adjacent subsamples were collected at each site to form a composite sample. Each sample campaign included two replicates and was stored at −20 • C for DNA extraction and molecular analyses.

Analyses of Physicochemical Parameters of Water
At each sample site, an automated multiparameter profiler (model YSI EXO) was used to measure water temperature (T), conductivity (Cond), total dissolved solids (TDS), pH, dissolved oxygen (DO), and chlorophyll a (Chla) in situ. Total nitrogen (TN) was tested using ion chromatography (EPA 300.0, revision 2.1, 1993) with prior persulfate oxidation for the water samples (Hosomi and Sudo, 1986). Samples for nitrate (NO 3 − ) and ammonium (NH 4 + ) analyses were acidified by addition of H 2 SO 4 (pH < 2) and stored at <4 • C in the dark prior to analysis. Concentrations of NO 3 − and NH 4 + were measured with an automatic flow analyzer (SKALAR Sans Plus Systems, Breda, Netherlands). NO 3 − was tested using ion chromatography (EPA 300.0, revision 2.1, 1993). The indophenol colorimetric method was used to analyze NH 4 + (EPA 350.1, revision 2.0, 1993). The ascorbic acid colorimetric method with prior oxidation used to determine total phosphorus (TP) (EPA 365.3, 1978). The main physical properties of each water sample are shown in Supplementary Table 1.

DNA Extraction and 16S rRNA Sequencing
Genomic DNA was extracted from the filter and sediment (0.5 g) samples using the FastDNA SPIN Kit for Soil (MP, United States), following the manufacturer's instructions. The V3 − V4 region of bacterial 16S rRNA genes was amplified using the forward primer 338F (5 -ACTCCTACGGGAGGCAGCAG-3 ) and reverse primer 806R (5 -GGACTACHVGGGTWTCTA AT-3 ) (Xu et al., 2016). The PCR was performed using the following protocol: initial denaturation at 94 • C for 5 min, followed by 25 cycles of denaturation at 94 • C for 30 s, annealing at 55 • C for 30 s, and extension at 72 • C for 30 s; and a final extension step at 72 • C for 5 min. Each sample was amplified in triplicate, pooled, and purified using the Agencourt AMPure XP Beads (Backman, United States). DNA amplicons were quantified using the Qubit TM dsDNA Assay Kit (Invitrogen, Carlsbad, CA, United States). Subsequently, an equi molar concentration of each DNA amplicon was mixed and then sequenced using the MiSeq PE300 platform (Illumina, San Diego, CA, United States) at Majorbio Bio-Pharm Technology Co. Ltd., Shanghai, China 1 .

Quantitative Real-Time PCR (qPCR)
The absolute abundance of bacterial 16S rRNA, ammonia monooxygenase of ammonia oxidizing bacteria (AOB), narG from nitrate reductase (narG), nirS from nitrite reductase (nirS), norB from reductase (norB), and nosZ from nitrous oxide reductase (nosZ) were measured by qPCR using bacteria-specific primer pairs. Primers used for the amplification of these genes are listed in Supplementary Table 2. A 10-fold serial dilution of standard DNA plasmids was prepared by cloning the target gene fragment amplified from samples using the same primer used for qPCR. The PCR reaction mixtures contained 10 µL of 2 × SYBR R Premix Ex Taq TM II (Takara, Japan), 1 µL (10 pmol/µL) of both forward and reverse primers, and 1 µL template DNA in a final volume of 20 µL. All qPCR procedures were performed on a CFX96 Real-Time PCR Detection System (Bio-Rad Laboratories, Hercules, CA, United States) with an initial denaturation at 95 • C for 2 min, followed by 40 cycles lasting 20 s each at 95 • C, 60 s at 63 • C, and 60 s at 68 • C. The specificity of each amplicon was checked by melt curve analysis when PCR amplification was completed. A standard curve was constructed by plotting the relative fluorescent units at a threshold fluorescence value (CT) versus the logarithm of copy number of the standard plasmid DNA. The efficiencies of qPCR ranged from 85 to 98% (R 2 > 0.997). Each sample DNA and standard DNA were analyzed in triplicate.

Data Analysis
For sequencing data, the forward and reverse reads were joined, assigned to samples, and truncated by splicing the barcode and primer sequence based on the overlap relationship between paired-end reads. The quality of reads and the effect of merging were filtered by QIIME (version 1.9.1 2 ) to obtain the optimized data. The effective sequences were grouped into operational taxonomic units (OTUs) at a 97% sequence identity level using Usearch (version 7.0 3 ). Subsequently, Mothur (version 1.30.2 4 ) was used to determine the richness index (Chao1 and ACE), coverage index (Good's coverage), phylogenetic diversity index (Faith's PD), and diversity index (Shannon and Simpson) for alpha diversity analysis. Taxonomic classification was assigned using the SILVA132 database (Release 132 5 ) to calculate the community composition of each sample. Based on the unweighted-UniFrac distance matrices, Principal coordinate analysis (PCoA) was created by Program R (version 3.0.1) to study the similarity or difference among samples. Analysis of similarity (ANOSIM) was performed to verify the significance of community compositions from different environmental categories of samples. The longest gradient lengths were less than 3.5 in preliminary detrended correspondence analysis. Thus redundancy analysis (RDA) was used to explore the relationship between bacterial community and environmental variables. The environmental factors subset was obtained by the maximum Pearson correlation coefficient of the distribution difference between the environmental factors and the sample communities. The samples and environmental factors were reflected on the same two-dimensional ranking chart by CANOCO software (version 4.5, Wageningen, Netherlands), and a Monte Carlo permutation test with 999 permutations was used to analyze the significance of the parameters (Zhang et al., 2021b). Based on bacterial 16S rRNA gene sequences, the metabolic functions were predicted using Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) version 1.1.0 6 . A co-occurrence network of bacterial community at the genus level was generated based on Spearman's correlation significant analysis (Zhang et al., 2021a,b). Calculation of the correlation between species established a connection between species at a certain threshold correlation coefficient. The species correlation network diagram was generated. One-way ANOVA was performed using SPSS Statistics 20 (IBM, United States). Graphs were generated using GraphPad Prism 5 software.

Accession Numbers
The 16S rRNA gene sequences from this study have been deposited in the public NCBI database 7 under the accession numbers PRJNA680428.

Vertical Variation of Physical and Chemical Water Indices
The water Chla concentration of HJD Reservoir ranged from 0.21 µg/L to 2.36 µg/L and corresponded with the characteristics of an oligotrophic reservoir. The water temperature decreased from 18.3 • C at the surface to 12.1 • C at the bottom (Supplementary Table 1). TN and TP concentrations 6 http://picrust.github.io/picrust/ 7 http://www.ncbi.nlm.nih.gov/ fluctuated from 3.74 mg/L to 4.05 mg/L and from 0.041 mg/L to 0.085 mg/L, respectively. WJD Reservoir was in a eutrophic state, with a maximum Chl-a concentration of 15.5 µg/L on the surface. The water temperature ranged from 12.8 to 20.2 • C. TN and TP concentrations followed similar trends to those of HJD Reservoir, ranging from 3.63 mg/L to 4.41 mg/L and from 0.065 mg/L to 0.296 mg/L, respectively. The DO concentration remained high throughout the water column in two reservoirs, varying from 6.35 mg/L to 9.42 mg/L and from 6.34 mg/L to 13.1 mg/L, respectively. The water bodies in reservoirs were stratified based on the measurement of water temperature (T) at different water depths. In HJD Reservoir, 0 m was the epilimnion, 5-15 m was the thermocline, and below 15 m was the hypolimnion. In WJD Reservoir, 0 m was the epilimnion, 5-15 m was the thermocline, and below 15 m was the hypolimnion.

Abundance of Bacteria
With increasing water depth, copy numbers of bacterial 16S rRNA genes slightly increased at 5 m and then slightly decreased at 30 m in HJD Reservoir, while the decrease was greater in WJD Reservoir (Figure 2). In reservoir profile water, the average 16S rRNA gene copy numbers varied from 2.38 × 10 7 to 1.74 × 10 9 copies/L. The copy numbers of bacterial 16S rRNA genes were 7.84 × 10 11 to 9.67 × 10 12 copies/L on average in sediments, which were based on the density of the solid at approximately 2.5 g/cm 3 (Xu et al., 2006) (Supplementary Figure 1). The lowest number of 16S rRNA gene copies was found in the hypolimnion, followed by the epilimnion and thermocline. The highest number was present in the sediment. The vertical distribution patterns of bacterial abundance in the two reservoirs were consistent with the Chla concentrations. The abundance of bacterial 16S rRNA genes in the sediment was significantly different from water (p < 0.05). The bacterial 16S rRNA copies were negatively correlated with D (p < 0.01), TN (p < 0.01), and NO 3 − -N (p < 0.01), and positively correlated with T (p < 0.01), DO (p < 0.01), pH-w (p < 0.01), and Chla (p < 0.01) in water samples ( Table 1).

Bacterial Diversity and Richness
Generally, 1,422,843 raw sequences were obtained after quality filtering for 26 libraries. The high quality reads of each sample ranged from 30,121 to 44,687 (Supplementary Table 3). Based on the 97% identity threshold (Figure 2), the OTU numbers in HJD and WJD reservoirs ranged from 317 to 536, and 369 to 884, respectively. The Chao 1 index and Faith's PD of water samples were 473-1,920, and 48.0-134, respectively. Both OTUs and the Chao 1 index were significantly higher in the HJD Reservoir than those in the WJD Reservoir. In water samples, the results in the Chao 1 index corresponded well with those of OTUs, and continued to increase along with the water depth of the reservoirs. In the epilimnion and thermocline of the HJD Reservoir, OTUs and Chao1 showed no significant change, but the hypolimnion exhibited significantly higher changes than the epilimnion and thermocline. However, in the WJD Reservoir, OTUs and Chao1 were not affected by the thermocline, but were inverse to the trend observed for bacterial abundance. The OTU numbers and Chao 1 index of sediments in the two reservoirs varied from 2938 to 3027, and 4,833 to 4,982, respectively, which were significantly higher (3-10 times) than those in the water samples. These measurements suggested that sediments harbored more diverse bacteria than the water within each site. Both water and sediment of WJD Reservoir had greater bacterial diversity and richness. Correlation analysis demonstrated that OTUs and Chao 1 index were positively correlated with D (p < 0.01), TN (p < 0.01), NO 3 − -N (p < 0.01), TP (p < 0.01), and PO 4 3− -P (p < 0.01) and negatively correlated with T (p < 0.01), DO (p < 0.01), and pH-w (p < 0.01) in water samples from the two reservoirs (Table 1).

Bacterial Community Structures and Environmental Variables
The bacterial community among the 26 samples was analyzed by PCoA based on the unweighted-unifrac matrix at the OTU level ( Figure 3A). The variance of 37.5% was attributed to PC1, whereas PC2 captured 13.4% of the variance. The whole water communities were found to be separated into six distinct groups and the sediment samples from the two reservoirs clustered together. Corresponding to the water stratification, the epilimnion, thermocline, and hypolimnion of the two reservoirs clustered together. Consistent with the results of bacterial abundance and diversity, significant differences in bacterial communities were observed between sediments and water as well as in different water layers.
The relative abundance of the bacterial community at the phylum level is shown in Supplementary Figure 3A. Overall, the two reservoirs contained 59 bacterial phyla, including 11 and 17 dominant phyla (relative abundance > 1% of the total in each sample) in the water column and sediment ( Supplementary  Figure 2), respectively. Proteobacteria, Actinobacteria, and Bacteroidetes were the predominant phyla across the water columns of the two reservoirs. In HJD Reservoir, Proteobacteria ranged from 22.31 to 78.04% and exhibited significant vertical differences. The relative abundance of Proteobacteria in bottom water had a high degree of outbreak during the stratification period, accounting for more than 78% of the bacterial community at 50 m. This value was much higher than that the values in the middle and upper waters, with reduced proportions of other phyla, including Actinobacteria, Bacteroidea, Cyanobacteria, and Verrucomicrobia. Actinobacteria represented 32.99-35.47% of the phyla, with no vertical difference in 0 to 30 m, but minimum proportion (9.60%) in 50 m. The relative abundance of Bacteroidetes decreased with increasing water depth. Chloroflexi, Acidobacteria, Nitrospirae, and Chlorobi were higher from 30 to 50 m than from 0 to 15 m. In WJD Reservoir, Actinobacteria, Bacteroidetes, and Verrucomicrobia represented 22.39 to 38.60%, 11.62 to 22.72%, and 1.67 to 4.21%, of the phyla, respectively, with no vertical difference. Cyanobacteria ranged from 1.11 to 16.46% with the maximum abundance (16.46%) at 0 m. The proportion of Planctomycetes, Chloroflexi, Armatimonadetes, Acidobacteria, Nitrospirae, and Chlorobi showed significant vertical differences and was higher from 15 to 60 m than from 0 to 5 m.
The bacterial classes were also investigated, and twenty dominant classes (relative abundance > 1%) were screened among samples (Figure 3B). In HJD Reservoir, both γand β-Proteobacteria exhibited significant vertical differences. γ-Proteobacteria increased with increasing water depth, with the maximum abundance (65.24%) at 50 m and minimum proportion (0.53%) at 0 m. β-Proteobacteria were more abundant at 5 and 15 m than at 0, 30, and 50 m. Sphingobacteriia, Flavobacteriia, Cytophagia, Armatimonadia, and Verrucomicrobia were higher from 0 to 5 m than from 15 to 50 m, while that of Planctomycetacia, Nitrospira, and Chlorobia was higher from 30 to 50 m than from 0 to 15 m. In WJD Reservoir, γand β-Proteobacteria showed inverse trends, accounting for 3.41-27.64% and 8.95-25.19%, respectively. δ-Proteobacteria were higher from 15 to 60 m than from 0 to 5 m and α-Proteobacteria showed no vertical difference. Flavobacteriia were higher from 0 to 5 m than from 15 to 60 m, whereas Anaerolineae, Nitrospira, and Acidobacteria were higher from 30 to 60 m than from 0 to 15 m. Phycisphacrae and Chlorobia abundance was higher from 15 to 60 m than from 0 to 5 m.
At the genus level, the dominant genera were defined as a relative abundance of more than 5% (Supplementary Figure 3B). In HJD Reservoir, the hgcl_clade was the dominant genus at 0, 5, and 15 m, hgcl_clade, Pseudomonas and CL500-29-marine-group were dominant at 30 m. The proportion of Pseudomonas increased with increasing water depth, especially at 50 m, and Pseudomonas was the absolute dominant genus, resulting in the proportional compression of the other genera. Limnohabitans showed no vertical differences. The abundances of Unclassified_f_Sporichthyaceae, Flavobacterium. and Polynucleobacter were higher from 0 to 5 m than from 15 to 60 m, whereas those of Norank_f_Anaerolineaccae were higher from 30 FIGURE 2 | Variations in bacterial abundance, operational taxonomic units (OTUs), Chao1, and water temperature (T) at different water depths in Hongjiadu (HJD) and Wujiangdu (WJD) reservoirs. Error bars represent standard deviation of the mean (n ≥ 2). Different letters above the columns indicate significant differences at p < 0.05.  Bolded values represent significant correlations between the variables (p < 0.05). *p < 0.05; **p < 0.01.
to 50 m than from 0 to 15 m. Unclassified_f_Comamonadaceae and Terrimicrobium were more abundant at 15 m than at other water depths, accounting for 14.05 and 6.45%, respectively. In WJD Reservoir, hgcl_clade, Limnohabitans, and Cyanobacteria were dominant at 0 m, hgcl_clade was dominant at 5 and 30 m, Acinetobacter was dominant at 15 and 45 m, and hgcl_clade and Limnohabitans were dominant at 60 m. Acinetobacter and Limnohabitans represented 2.86-25.08% and 2.22-17.08%, respectively, with significant vertical differences. Unclassified_f_Sporichthyaceae, Flavobacterium and unclassified_f_Comamonadaceae were more abundance from 0 to 5 m than from 15 to 60 m. The distributions of the bacterial community at the phylum, class, and genus levels clearly exhibited significant vertical differences in the stratification of reservoir water.

Relationship Between Bacterial Community and Environmental Variables
Redundancy analysis (RDA) evaluated the influences of environmental factors on bacterial communities in water ( Figure 3C). The eigenvalues of the first two principal components closely approached 73.7%. The cumulative contribution ratios of RDA 1 and RDA 2 were 47.5 and 26.2%, respectively. The results indicated that TN (p < 0.01), NO 2 − -N (p < 0.01), NO 3 − -N (p < 0.05), TP (p < 0.05), and PO 4 3− -P (p < 0.05) were significantly positively correlated with the variations in bacterial communities in water. Overall, nutrient nitrogen was a crucial driver in the water bacterial communities. Based on the Pearson correlation analysis, the relationship between environmental parameters and the bacterial community was investigated (Figure 3D and Supplementary Figure 4). At the phylum level, Bacteroidetes, Cyanobacteria, Chloroflexi, Armatimonadetes and Nitrospirae were significantly correlated with T. Actinobacteria, Armatimonadetes, TM6_Dependentiae_, Nitrospirae, Chlorobi, and Acidobacteria were significantly correlated with DO (Supplementary Figure 4A). At the class level, Flavobacteriia, Phycisphaerae, Anaerolineae, and Cytophagia were significantly positively correlated for NO 3 − -N, NH 4 + -N, TN, PO 4 3− -P, and TP. γ-Proteobacteria and β-Proteobacteria exhibited a negative relationship with  Figure 4B). At the genus level, Limnohabitans, Fluviicola, Acinetobacter, and Sediminibacterium showed a positive relationship with age, but a negative relationship with HRT ( Figure 3D). Pseudomonas and Aquabacterium were significantly positively related with HRT, but negatively related with age. Hgcl_clade, norank_c_Cyanobacteria, unclassified_f_Sporichthyaceae, Limnohabitans, Polynucleobacter, Flavobacterium and Fluviicola showed a positive relationship with T.

Prediction of Microbial Metabolism and Co-occurrence Network Analysis
To understand the biological features and functional contributions of the microbial population at different water depths, metagenomes were predicted by PICRUSt using the 16S rRNA gene amplicon data. Metabolism pathways of microbial community predicted by PICRUSt were further analyzed as shown in Figure 4. The relative abundance of the metabolism pathway was the highest among the level 1 Kyoto Encyclopedia of Genes and Genomes (KEGG) categories in two reservoirs, accounting for 73.0-78.9% in HJD Reservoir and from 77.4 to 78.4% in WJD Reservoir (Figure 4A). Carbohydrate metabolism and amino acid metabolism were the top two pathways among the metabolism categories at level 2 ( Figure 4B). Carbohydrate metabolism including pyruvate metabolism, glyoxylate and dicarboxylate metabolism, propanoate metabolism, butanoate metabolism, and others, was relatively highly abundant at 50 m in HJD Reservoir. At the same depth, amino acid metabolism pathways were also highly abundant. These included glycine, serine and threonine metabolism, cysteine and methionine metabolism, valine, leucine and isoleucine degradation, alanine, aspartate and glutamate metabolism, arginine and proline metabolism, and others. The top 30 metabolic pathways among the level 3 KEGG categories are listed in Figure 4C. Generally, the biosynthesis of amino acids and carbon metabolism were the most abundant pathways in the two reservoirs. Interestingly, environmental information processing related pathways, such as ABC transporters and two-component system, displayed high relative abundance in HJD Reservoir compared to that in WJD Reservoir, especially at 50 m, indicating different enriched KEGG pathway patterns between these reservoirs.
Co-occurrence network analysis clarified the coexistence relationship among different bacterial genera ( Figure 5A). Pseudomonas and Aquabacterium exhibited a significant symbiotic relationship. These two genera were significantly repellent to other dominant bacterial genera. In addition, a significant symbiotic relationship was evident between unclassified_Cyanobacteria, hgcI_clade, unclassified_Comamonadaceae, Sediminibacterium, Limnoh abita, Fluviicola, Polynucleobacter, unclassified_Sporichthyaceae and Flavobacterium. Correlations between the bacterial relative abundance and the enriched KEGG orthologous (KO) gene families are presented in Figure 5B. In brief, 41 KOs were significantly positively correlated with Pseudomonas abundance (p < 0.01 or 0.05), including K02003 (putative ABC transport system ATP-binding protein) and K02004 (putative ABC transport system permease protein) involved in the environmental information processing related pathways. Unclassified_Cyanobacteria and hgcI_clade showed positive correlations with K01952 (phosphoribosylformylglycinamidine synthase), K03657 (DNA helicase II/ATP-dependent DNA helicase PcrA) and K01358 (ATP-dependent Clp protease, protease subunit) (p < 0.01 or 0.05). The abundance of Sediminibacterium and Fluviicola was negatively associated with 23 KOs and 11 KOs, respectively (p < 0.01 or 0.05). Fluviicola, unclassified_Sporichthyaceae and Flavobacterium exhibited positive correlations with both K01952 (phosphoribosylformylglycinamidine synthase) and K00826 (branched-chain amino acid aminotransferase). CL500-3 and CL500-29 marine group exhibited strong negative correlations with K01996 (branched-chain amino acid transport system ATP-binding protein) and K01998 (branched-chain amino acid transport system permease protein), respectively (p < 0.01 or 0.05).

Absolute Abundance of Genes Associated With Major Nitrogen Metabolism Pathways
In HJD Reservoir, the gene abundance of the nitrification gene (AOB) ranged from 9,731 to 1.02 × 10 5 copies/L ( Figure 6A). The AOB abundance decreased from 0 to 15 m, and then increased from 15 to 50 m. However, these changes were not statistically significant. Gene copies of narG and nosZ genes ranged from 9.50 × 10 4 to 3.21 × 10 7 copies/L and from 1.55 × 10 4 to 5.99 × 10 6 copies/L, respectively. The abundance of narG and nosZ genes increased from 0 to 15 m, decreased sharply at 30 m, and then significantly increased at 50 m. The abundance of nirS and norB genes ranged from 2.14 × 10 5 to 6.88 × 10 6 copies/L and from 3.31 × 10 4 to 8.24 × 10 5 copies/L, respectively. No significant difference in the nirS and norB genes was observed among the water layers. In WJD Reservoir, the copy number of the AOB gene ranged from 1.16 × 10 4 to 3.83 × 10 6 copies/L. The temporal changes in the abundance of this gene were similar to those observed in HJD Reservoir. The copy numbers of narG and nosZ genes ranged from 5.44 × 10 5 to 1.76 × 10 7 copies/L and from 1.03 × 10 4 to 3.07 × 10 6 copies/L. The abundance was greatest at 0 m and then decreased from 5 to 60 m (although not significantly). The copy numbers of nirS and norB genes ranged from 1.40 × 10 5 to 2.97 × 10 7 copies/L and from 8.80 × 10 4 to 1.18 × 10 7 copies/L, respectively. The abundance of nirS and norB genes was the highest at 0 m, and no significant difference was observed between 5 and 60 m. Correlation analysis (Figure 6B) revealed the key physicochemical factors explaining changes in the expression of all selected genes. AOB abundance was significantly positively correlated with PO 4 3− -P and TP; narG and nosZ copy numbers were significantly negatively correlated with NH 4 + -N and positively correlated with NO 2 − -N; nirS and norB gene copies were significantly negatively correlated with NO 3 − -N and TN.

Spatial Dynamic Distribution of Bacterial Abundance in Reservoirs of Different Trophic Statuses
The trends in variation of bacterial abundance did not correspond with different trophic reservoirs. The spatial change in bacterial abundance in HJD Reservoir was relatively slight. HJD Reservoir was reported as an oligotrophic reservoir (Zhao et al., 2017) with long hydraulic residence time and high light intensity. Bacteria selectively grew in the water layer with appropriate light, and displayed increased abundance with increasing water depth in the epilimnion layer. The appearance of a thermocline reduces the levels of particulate matter (Wilhelm and Adrian, 2008), which provides an epiphytic habitat for microorganisms, resulting in bacterial accumulation. The lack of nutrients and decrease in DO vertical migration (Yu et al., 2014;Zhang et al., 2015) limited the bacterial growth in the hypolimnion layer. A remarkable spatial dynamic distribution of bacterial abundance has been observed in WJD Reservoir (Figure 2). WJD Reservoir has been facing a serious threat to eutrophication (Zhao et al., 2017). The trophic status might be a key driving force for increased bacterio-plankton abundance in freshwater (Lindstrm, 2000;Dai et al., 2016). Excessive nutrients and weak light transmittance can promote the growth and reproduction of phytoplankton in surface water (Rother, 1977;Feng et al., 2019). In this study, Chla and cyanobacteria contents were 15.5 µg/L and 16.5%, respectively, WJD surface water. Furthermore, the FIGURE 6 | Abundance of nitrogen functional genes and the correlation of environmental factors. (A) The absolute abundance of nitrification and denitrification genes in water stratification of the two reservoirs. Error bars represent the standard deviation (n ≥ 4). Different letters indicate significant differences at p < 0.05 level (one-way ANOVA). (B) SPSS correlation analyses of the abundance, diversity, and environmental nutrient factors. Arrows represent positive (red) or negative (green) correlation between the variables, respectively. The arrow width is proportional to the correlation strength. Numbers on the arrows are the correlation coefficients. * p < 0.05; * * p < 0.01. growth and decomposition of phytoplankton could induce the growth of bacteria (Paver et al., 2013). This may be the reason for the higher bacterial abundance in surface water of the more eutrophic WJD Reservoir than that in the surface water of the oligotrophic HJD Reservoir. However, excessive oxygen consumption by algae and the release of harmful algal toxins (Berg et al., 2009) might lead to extensive bacterial death (Nitin et al., 2017). In addition, the thermocline hindered the mixing of the water column, limited the vertical migration of DO and particulate matter, increased the heterogeneity of the water column (Selmeczy et al., 2016), and resulted in the generation of an environment that was less conducive to bacterial growth. Meanwhile, low DO stimulated the release of nutrients from the sediment to the overlying water (Rigosi and Rueda, 2012), which would worsen the water quality (Elçi, 2008) and might aggravate the degree of eutrophication of WJD Reservoir.

Thermal Stratification and Nutrients Might Play Key Roles in Shaping the Vertical Dynamic Distribution of the Bacterial Community
In this study, OTUs, richness, and diversity of bacteria in the sediments of the two reservoirs were found to be significantly higher than those in the water, consistent with the results of previous studies (Dai et al., 2016;Feng et al., 2019;Yue et al., 2021). In the vertical distribution of the water column, the OTUs, richness, and diversity of bacteria in eutrophic reservoirs were also significantly higher than those in oligotrophic reservoirs and were significantly correlated with variables related to the nutritional status. Recent studies have shown that the nutritional status might be an important driving force for increased bacterial community diversity (Li et al., 2017;Sun et al., 2020). In our study, bacterial OTUs, richness, and diversity of the hypolimnion were found to be significantly greater than those of the epilimnion and thermocline in the same reservoir. It was possible that the thermal stratification caused hypoxia in the hypolimnion and stimulated nutrient decomposition in the sediment, which increased the number of microorganisms entering the hypolimnion.
Actinobacteria (8.77-37.3%) was the dominant phylum at most water depths in the two reservoirs (Allgaier and Grossart, 2006). The finding indicates that Actinobacteria was unaffected by the nutritional status and mixed water environment (Nitin et al., 2017). However, compared to WJD Reservoir (the eutrophic reservoir), HJD Reservoir (the oligotrophic reservoir) displayed a more pronounced vertical spatial difference in terms of bacterial community structure during thermal stratification. Correspondingly, predictive metagenomic profiling revealed different vertical functional variation patterns of the bacterial community in the two reservoirs. Bacteroidetes (26.1-28.6%), β-Proteobacteria (16.7-28.1%), and γ-Proteobacteria (16.8-66.2%) were dominant in the epilimnion, thermocline, and hypolimnion of HJD Reservoir, respectively. Pseudomonas was the most dominant genus at 50 m, and was significantly involved in 44 KOs ( Figure 5B and Supplementary Table 4), such as the functional genes related to metabolism and environmental information processing pathways. Considering that the competition among the bacterial communities for the nutrient resource might be intense, given the oligotrophic nature of this reservoir, the high abundance of Pseudomonas at 50 m might be associated with the multiple roles of these bacteria in the water column.
Nutritional status is considered an important determinant of various aspects of bacterial community structure (Wu et al., 2010;Jankowski et al., 2014;Davis et al., 2015;Zwirglmaier et al., 2015). It has been reported that TP and TN significantly contribute to the bacterial community in the Qingliu River that had been affected by domestic sewage (Zhang et al., 2020). Li et al. (2017) also showed that the nutrients (especially total carbon) play a more important role rather than the sample locations in shifting bacterial structure in the shallow small lakes. Our RDA results showed that the bacterial community structure in the two reservoirs was significantly correlated with nutrient parameters, in agreement with the results from previous studies. Furthermore, the co-occurrence network analysis based on our dataset revealed the significant symbolic relationships among Cyanobacteria, Actinobacteria (hgcl_clade, Sporichthyaceae) and β-Proteobacteria (Comamonadaceae, Limnohabita, Polynucleobacter), consistent with the reported observations from other studies. For instance, α-Proteobacteria (Berg et al., 2009), β-Proteobacteria (Eiler et al., 2012) and Actinobacteria (Zhao et al., 2016) are respectively in physical association with some members of Cyanobacteria in freshwater systems. Similarly, Chloroflexi, Acidobacteria, Cyanobacteria, and Planctomycetes have significant co-occurrence relationships with many phyla in the surface water in Xukou Bay of Lake Taihu (Cao et al., 2018).

Vertical Dynamic Distribution of Nitrification and Denitrification in the Water Column
In this study, qPCR was used to estimate the vertical distribution of water column nitrification and denitrification abundance in HJD and WJD reservoirs. In HJD Reservoir, nitrification and denitrification in the hypolimnion were much higher than those in the thermocline. DO was the main limiting factor in denitrification owing to the lack of nutrients in the water column. The emergence of a thermocline reduced the DO content and increased denitrification. This occurred because the DO content in the upper water could not supply the hypolimnion, which stimulated an increase in denitrification. Despite the general assumption that aerobic conditions are conducive to nitrification, nitrification was greatest in the hypolimnion. Correlation analysis of the process of nitrification revealed a positive correlation between amoA and nirS and norB. These findings suggested that nitrification and denitrification had coupled and exhibited synergistic effects in the water column. In the hypolimnion, denitrification resulted in the production of inorganic substances and provided NH 4 + to facilitate nitrification. Meanwhile, with the increased nitrification in the hypolimnion, DO was continuously consumed and promoted denitrification.
However, in WJD Reservoir, no significant difference was observed in the nitrification and denitrification processes in the different thermal strata of the water column, indicating that the sensitivity of nitrification and denitrification was less affected by heterogeneity in the vertical direction of the reservoir's water column. The correlation analysis also revealed that NO 3 − -N, NH 4 + -N, NO 2 − -N, TN, TP, and PO 4 3− -P were significantly correlated with nitrification gene (AOB) and denitrification genes (nirS, norB). These findings suggest that nutrient sufficiency is the main limiting factor for nitrification and denitrification processes in eutrophication reservoirs (Tang et al., 2016). The sudden reduction in DO (Yu et al., 2014;Zhang et al., 2015) and the slow migration of particulate matter (Wilhelm and Adrian, 2008) caused by thermal stratification had little effect on nitrification and denitrification.

CONCLUSION
We examined the bacterial community structures based on quantitative real-time PCR and high-throughput sequencing in two different trophic reservoirs. During water stratification in the water column, bacterial abundance in HJD Reservoir was slightly affected by the thermocline. In WJD Reservoir, bacterial gene copies sharply decreased with increasing water depth and exhibited obvious spatial differences. Proteobacteria, Actinobacteria, and Bacteroidetes were the dominant phyla. The dominant classes and genera were influenced by changes in water stratification and trophic levels. RDA and heatmap analysis indicated that thermal stratification and nutrients were derived from variations in the bacterial community structure. Moreover, significant differences were evident in nitrification and denitrification. The latter was markedly more extensive than nitrification, regardless of nutritional status. These results provide valuable insights into mechanisms by which water stratification and trophic levels control bacterial abundance and structure in cascade reservoirs.
A limitation of the study was that it was conducted in April and stratification is more pronounced during summer. A study conducted in summer may reveal more interesting variations in bacterial structure.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The name of the repository and accession number can be found below: National Center for Biotechnology Information (NCBI), https://www.ncbi.nlm.nih. gov/, PRJNA680428.

AUTHOR CONTRIBUTIONS
All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.