Shifts in the Microbial Community, Nitrifiers and Denitrifiers in the

The objective of this study was to investigate the microbial community shifts, especially nitrifiers and denitrifiers, in the biofilm of two rotating ...
0 downloads 0 Views 2MB Size
Download PDF
Article pubs.acs.org/est

Shifts in the Microbial Community, Nitrifiers and Denitrifiers in the Biofilm in a Full-scale Rotating Biological Contactor Xingxing Peng, Feng Guo, Feng Ju, and Tong Zhang* Environmental Biotechnology Laboratory, Department of Civil Engineering, The University of Hong Kong SAR, China S Supporting Information *

ABSTRACT: The objective of this study was to investigate the microbial community shifts, especially nitrifiers and denitrifiers, in the biofilm of two rotating biological contactor (RBC) trains with different running times along the plug flowpath. The microbial consortia were profiled using multiple approaches, including 454 high-throughput sequencing of the V3−V4 region of 16S rRNA gene, clone libraries, and quantitative polymerase chain reaction (qPCR). The results demonstrated that (1) the overall microbial community at different locations had distinct patterns, that is, there were similar microbial communities at the beginnings of the two RBC trains and completely different populations at the ends of the two RBC trains; (2) nitrifiers, including ammonia-oxidizing archaea (AOA), ammonia-oxidizing bacteria (AOB, Nitrosomonas) and nitrite-oxidizing bacteria (NOB, Nitrospira), increased in relative abundance in the biofilm along the flowpath, whereas denitrifiers (Rhodanobacter, Paracoccus, Thauera, and Azoarcus) markedly decreased; (3) the AOA were subdominant to the AOB in all sampled sections; and (4) strong ecological associations were shown among different bacteria. Overall, the results of this study provided more comprehensive information regarding the biofilm community composition and assemblies in full-scale RBCs.

■

INTRODUCTION Nitrogen removal, including the processes of nitrification and denitrification is one of the most important targets of wastewater treatment plants. The rotating biological contactor (RBC) is a system that uses immobilized microorganisms in biofilm as catalysts to trigger a series of biological reactions to remove pollutants, including nitrogen.1−6 These systems have been widely used because of their reduced land demand and susceptibility to temperatures and their lower energy consumption compared with other treatment processes. In the RBC, biofilms attached to the surface of discs and the assimilation of the pollutants in the wastewater play the dominant roles during wastewater treatment. Recent studies have reported several specific species of nitrifiers and denitrifiers;7−12 however, the biofilm community composition in RBCs has been poorly characterized compared with other biological wastewater treatment systems. In recent decades, several molecular technologies based on microbial 16S rRNA, including fluorescence in situ hybridization (FISH),7−10 denaturing gradient gel electrophoresis (DGGE),9,11 qPCR12 and clone libraries have been used to profile the nitrogen-removing bacteria and organics-removing bacteria in RBC biofilms.8 Sauder et al.12 used quantitative PCR and PCR-DGGE technology to investigate ammonia-oxidizing populations, including ammonia-oxidizing bacteria (AOB) and archaea (AOA) in RBC biofilms from a municipal wastewater treatment plant, and observed that an ammonia gradient determined the relative abundance of AOA and AOB. Recently, © 2014 American Chemical Society

454 pyrosequencing, generating large numbers of intermediatelength DNA reads through a massively parallel sequencing,13 has been used to comprehensively determine the bacterial community in different environmental samples, including ocean, soils, and activated sludge.14−16 In this study, biofilm samples were collected from different parts of a full-scale RBC at the Ma Wan Sewage Treatment Plant in Hong Kong. The pyrosequencing of 16S rRNA gene was conducted to analyze the bacterial biofilm community at the level of phylum, genus and species based on a linear sampling strategy to profile the variation in the bacterial consortia along the RBC flowpath. In addition, qPCR and clone library analyses were conducted to detect the relative abundance and diversity of AOA and AOB amoA sequences. Furthermore, the similarity/difference among different sections of the RBC was evaluated, and the correlations of those functional bacteria were identified. As the first application of pyrosequencing in this context, a more comprehensive analysis of the bacterial community in the RBC was unraveled, and a shifting pattern was found along with the pollutant removal. Received: Revised: Accepted: Published: 8044

April 7, 2014 June 17, 2014 June 17, 2014 June 17, 2014 dx.doi.org/10.1021/es5017087 | Environ. Sci. Technol. 2014, 48, 8044−8052

Environmental Science & Technology

Article

Figure 1. Schematic diagram of rotating biological contactor, sampling name and location. The samples RBC-A-1 to RBC-A-7 were collected from the second set (Set A), and the samples RBC-B-1 to RBC-B-7 were collected from the fourth set (Set B). Samples from RBC-A-1a to RBC-A-1c and from RBC-B-7a to RBC-B-7c were duplicates located in the 1st disc of Set A and the 7th disc of Set B, respectively.

■

MATERIALS AND METHODS Sample Collection. Eighteen samples were collected from two trains of rotating biological contactors (Set A and Set B) in Ma Wan Sewage Treatment Plant in Hong Kong, which treats, on average, 7600 m3 per day of municipal wastewater (Figure 1). Set A and Set B samples were collected from all stages of both the set A and set B treatment trains. Set A had been running for 6 months, whereas Set B had been running for 1 month. All biofilms samples were collected by aluminum metal panels (15 × 20 × 0.2 mm) and then immersed into 50 mL centrifuge tubes with 70% ethanol. RBC-associated water samples were also collected for each RBC stage and placed into sterile plastic tubes on ice (Figure 1). In addition, a mixture of all biofilm samples was prepared for the construction of clone libraries. Water Chemistry Analysis. TOC (total organic carbon) was determined using a Total Organic Carbon Analyzer (Shimadzu TOC-VCPH, Kyoto, Japan). COD (chemical oxygen demand) was measured by a standard potassium dichromate titration method.17 NH4+-N and NO3−-N were determined according to the Standard Methods for Water and Wastewater Examination. The pH value was determined using a pH meter (Thermo Orion 4 Star, Beverly, USA). DNA Extraction, PCR Amplification and Pyrosequencing. All biofilm samples were washed twice with a 0.85% NaCl solution at 4000 rpm for 12 min at 4 °C, and DNA was extracted using the FastDNA SPIN Kit for Soil (Q-Biogene, CA).18 The V3−V4 region (∼420 bp) of the bacterial 16S rRNA genes was amplified using the following primers: forward 5′- ACTCCTACGGGAGGCAGCAG-3′ and reverse 5′TACNVGGGTATCTAATCC-3′. The barcodes for multi-

plexing were inserted between the forward primers and the 454 adapter sequence.19 The reactions were run on an i-Cycler (BioRad Hercules, CA) under the following cycling conditions: a 5 min initial denaturation at 94 °C, followed by 30 cycles of denaturing at 94 °C for 30 s, annealing at 50 °C for 45 s, extension at 72 °C for 60 s, and a final extension at 72 °C for 10 min. After PCR amplification, the PCR products of three replicates were mixed and purified with a PCR quick-spin Kit (iNtRON Biotechnology, Korea). A NanoDrop device (ND1000) and electrophoresis were utilized to quantify the DNA concentration and to verify the fragment lengths, respectively. The PCR products were then sent to the Genome Research Centre at the University of Hong Kong for pyrosequencing conducted on a the Roche 454 FLX Titanium platform (Roche, Nutley, NJ). Processing of Pyrosequencing Data. The raw sequences from pyrosequencing were processed based on the Mothur v 1.25.0 analytic pipeline.20 After demultiplexing, quality trimming and alignment, chimeric sequences were removed after finally being checked with a ChimerSlayer algorithm. The remaining sequences were grouped into operational taxonomic units (OTUs) at genus and species levels using at 97% and 98% identity thresholds using the Quantitative Insights Into Microbial Ecology (QIIME) toolkit, respectively.21 These sequences were also assigned to corresponding taxonomic ranks using the RDP classifier. A bootstrap cutoff (80%) suggested by the RDP was applied to assign the sequences to different taxonomy levels. Good’s coverage was used to measure the captured diversity and was calculated as G = 1 − n/N, where n is the number of single OTUs, and N represents the total number of sequences. 8045

dx.doi.org/10.1021/es5017087 | Environ. Sci. Technol. 2014, 48, 8044−8052

Environmental Science & Technology

Article

Table 1. Summary of the Wastewater Parameters Obtained From the RBC Water Samplesa sample a b c d e a

TOC (mg/L) 80.2 49.6 35.2 28.6 15.9

± ± ± ± ±

5 3 1 2 1

NH4+-N (mg/L)

COD (mg/L) 327.4 202.9 130.5 100.4 47.7

± ± ± ± ±

10 7 5 6 2

115.2 93.8 69.7 37.2 4.6

X ± s: means ± standard deviation (n = 3).

± ± ± ± ±

5 4 2 4 1

NO3‑-N (mg/L) 2.5 1.8 5.6 8.7 17.7

± ± ± ± ±

1 1 1 2 2

pH 6.6 6.2 5.7 6.6 6.8

± ± ± ± ±

0.1 0.1 0.2 0.1 0.2

■

RESULTS AND DISCUSSION Pollutant and Physicochemical Gradients along RBCs. Throughout the whole stage, the concentrations of organic carbon and NH4+-N were generally decreased along the RBC flowpaths (Table 1). Overall, acceptable treatment capacities were found at the effluent of the RBC trains, with the removal of TOC (80%), COD (85%), and NH4+-N (96%). These results are consistent with previous studies that have shown the removal efficiency along the RBC flowpath.4−6 The pH value was observed to decrease and then increase. The decrease in pH value in RBC stage may be caused by the alkalinity consumption of the ammonium oxidizing process, whereas the increase in pH in the later RBC stage should be due to the alkalinity production.32 The denitrification process may reduce in the later RBC stage with the increase of NO3−-N (5.6−17.7 mg/L). Overall Community Analysis. After the cleaning of raw reads, a total of 49 918 effective sequences were obtained. To evaluate the validity of our test, the reproducibility of pyrosequencing was evaluated by computing the Pearson’s correlation coefficient. The results (Pearson’s correlation coefficients ranging from 0.900 to 0.971) show that the replicated samples at both the initial and end sections were technically highly reproducible.22,23 The sequencing depth was sufficient to describe patterns (Supporting Information (SI) Figure S1), and the bacterial coverage slightly decreased as the Good’s value declined from 98.3% to 97.7%. All of these results indicated that the library for each sample was fully satisfactory to characterize bacterial communities. These 49 918 effective bacterial sequences belonged to different taxa levels (from genus to phylum) with the RDP Classifier at an 80% threshold (SI Figure S2). A portion of the effective bacterial sequences could not be assigned to any taxa of different levels, which illustrates that some bacteria in the RBC biofilm were novel. Because the environment impacted bacterial composition,14 both the novel sequences and the bacterial diversity decreased along the RBC flowpaths as the pollutant concentrations gradually decreased. The relative abundances of different phyla and classes of Proteobacteria were analyzed (SI Figure S3), and in these 18 biofilm samples, Proteobacteria were the most abundant phylum, accounting for 34.8−57.6% of the total effective bacterial sequences. Within Proteobacteria, Alphaproteobacteria, Epsilonproteobacteria, Gammaproteobacteria, Betaproteobacteria, and Deltaproteobacteria were arranged in descending order of relative abundance. The other dominant phyla were Actinobacteria (3.9−24.5%, averaging 14.3%), Bacteroidetes (2.9−20.7%, averaging 12.0%), Nitrospirae (0.5−28.2%, averaging 5.4%), Firmicutes (0.9−6.7%, averaging 3.6%), and Chlorof lexi (0.7−5.6%, averaging 2.5%), which presented similar compositions at phylum level with other rotating biological contractors.33 The nitrite-oxidizing bacteria Nitrospirae were also detected and common along the RBC flowpath especially in the later RBC stage, such as the

The sequencing reproducibility of the pyrosequencing method was evaluated by assessing Pearson’s correlation coefficient22,23 using SPSS software. The target reproducibility of our study was calculated by finding a sufficient correlation coefficient of at least 0.8.24 Cloning and Sequencing of AOA and AOB amoA Genes. Clone libraries to investigate amoA diversity (including AOA and AOB) in the RBCs were constructed by providing longer sequence reads than those offered by 454 pyrosequencing. The primers were described in SI Table S1, and the experiments were conducted following a previous study.25 A total of 50 clones were sequenced for both AOA and AOB. Quantitative PCR (qPCR). The qPCR was performed in a 25 μL reaction that contained 12.5 μL of Maxima SYBR Green qPCR Master Mixes (Thermo), 10 pmol of each primer (SI Table S1), 20 ng of DNA template and ultrapure sterile water. Each sample had three replicates, and all qPCR was performed on an iCycler thermal cycler (Bio-Rad, Philadelphia, PA) equipped with iCycler iQ software (Bio-Rad). The PCR protocol was 1 min at 95 °C; 45 cycles consisting of 10 s at 95 °C, 20 s at 56 °C, and 1 min at 72 °C, then 50 s at 72 °C. The PCR conditions were the same for the AOB amoA genes except for the annealing temperature (58 °C). The results were effective based on the correlation coefficient of the standard curve (R2 = 0.995). Statistical Analysis. Cluster analysis (CA) and UniFrac principal coordinate analysis (PCoA) are statistical techniques that sort observations into similar groups or sets. CA was based on taxonomy results and operational taxonomic units (OTUs), whereas PCoA was based on RDP Classifier results, OTUs and UniFrac, which is a phylogeny-dependent method using phylogenetic information to compare different samples.26 Both of these analyses were conducted using PAST software.27,28 Network analysis of bacterial co-occurrence patterns may effectively reveal potential bacterial interactions in environmental samples.29 In this study, all pairwise Spearman’s rank correlations between 1949 OTUs (screening from 6457 OTUs) were calculated, and 370 OTUs were left after filtering, retaining only those with a Spearman’s correlation coefficient ρ > 0.6 and P-value < 0.01.30 This previous filtering step removed poorly represented OTUs as well as reduced the network complexity, facilitating the formation of the core community. The nodes represented the OTUs corresponding to 97% identity in the reconstructed network, whereas the edges depicted a strong and important correlation between nodes. Some other parameters, such as average node connectivity, diameter, average path length, clustering coefficient, and the modularity index, were also calculated. All statistical analyses were conducted using R with operating functional and positive packages, and the Gephi platform was used to explore and visualize networks.31 8046

dx.doi.org/10.1021/es5017087 | Environ. Sci. Technol. 2014, 48, 8044−8052

Environmental Science & Technology

Article

Figure 2. Heat map of genera (abundance >0.15%) in 18 biofilm samples. The color bar indicates the range of the percentage of a genus in a sample, based on the color key (log scale) at the top right. Genera with abundance >0.15% were selected in each sample (a total of 370 genera). Clustering of the cases was based on the genera abundance (>0.15%) data shown in the left heat map. Ammonia-oxidizing bacteria (AOB): red ★; Nitriteoxidizing bacteria (NOB): ▲; Denitrifier: ●.

sample (RBC-B-7), with the treatment at each stage creating a pollutant gradient along the flowpath. A heat map generated with major genera in 18 biofilm samples based on RDP analysis is shown in Figure 2. Genera with relative abundance >0.15% (48 genera) in each sample were selected and compared with relative abundances in other samples. The 48 genera belong to 7 phyla: Acidobacteria, Bacteroidetes, Chlorof lexi, Firmicutes, Proteobacteria, Spirochaetes, and TM7. Nitrosomonas (AOB) and Nitrospira (NOB) notablely increased along the RBC flowpath; it is likely that the water close to the effluent contained lower concentrations of pollutants (such as TOC, ammonia),37 consistent with the previous study.12 Four denitrifying groups, Rhodanobacter, Paracoccus, Thauera, and Azoarcus were also the dominant genera, which were more abundant in the later stages. In addition, the cluster analysis suggested that five clusters of genera were involved in this study (Figure 2, left). The five clusters represent stage-dependent changes in relative abundances of some of the dominant bacteria of in samples and experienced shifts along the wastewater treatment flowpath. A significant increase in the relative abundance of cluster I (Micromonospora, Arenibacter, Pelagibius, Rhodanobacter, Virgibacillus, Legionella, Nitrosomonas, and Erythrobacter) was detected from the RBC inlet to the outlet in both Set A and Set B. Cluster II (Nitrospira) also notablely increased in relative abundance, consistent with the enhancement of the denitrification process. In contrast, cluster III (Geminicoccus,

duplicates of RBC-B-7, which were notablely enriched along the RBC flowpath. These results are consistent with a previous study in that the structure and composition of the RBC biofilms correlated well with ammonium conversion.8 At the genus level, a total of 6457 OTUs in the 18 biofilm samples were constructed using 3% dissimilarity as the cutoff, and the top 20 OTUs were selected and numbered sequentially according to their relative abundance level from high to low (SI Figure S4). These top 20 OTUs belong to 6 phyla, including Proteobacteria (α-, β-, and γ-), Bacteroidetes, Actinobacteria, Nitrospirae, Chlorof lexi, and Firmicutes. The most abundant genus (including OTU1, OTU11, and OTU18), affiliated with Mycobacterium, has been reported to degrade many carboncontaining compounds, such as polycyclic aromatic hydrocarbon, phenanthrene, and pyrene, etc.34,35 From the heat map of the OTU relative abundances (SI Figure S4), OTU1, OTU11, and OTU18 all decreased along the RBC stages, and this trend was consistent with the decrease of TOC and COD along the RBC stages. The second highest genus, that is, Nitrospira (including OTU2 and OTU8), notablely increased along the RBC flowpath, especially in the later RBC stages. The shift in OTU 12 (Nitrosomonas, AOB), which was consistent with other studies,36 showed that AOB increased by 8.5 times (Set A) and 23 times (Set B) from the RBC inlet to the outlet, respectively. Overall, the total genera decreased from 169 in the initial biofilm sample (RBC-A-1) to 52 in the final biofilm 8047

dx.doi.org/10.1021/es5017087 | Environ. Sci. Technol. 2014, 48, 8044−8052

Environmental Science & Technology

Article

Figure 3. Variation in OTUs classified as Nitrosomonas and Nitrospira along the RBC flowpath after normalization by the most abundant OTU. The ratios were the relative abundance of OTUs divided by the OTU 1 (i.e., the most abundant OTU of Nitrosomonas or Nitrospira). The relative abundance for each specie was calculated by assigning sequences to the taxa and comparing them with the number of total 16S rRNA gene sequences in the RBC samples.

Figure 4. Phylogenetic tree of the archaeal amoA (a) and bacterial amoA (c) gene sequences retrieved from the biofilm samples along the RBC stages. Both trees were inferred using the maximum likelihood approach, and the bootstrap values are based on 1000 replicates. The scale bars represent 0.05 substitutions per site. The abundances of archaeal amoA (b) and bacterial amoA (d) were significantly increased along the RBC treatment process. 8048

dx.doi.org/10.1021/es5017087 | Environ. Sci. Technol. 2014, 48, 8044−8052

Environmental Science & Technology

Article

Paracoccus, Clostridium_XI, Mycobactenum, Rhodobacter, and Hydrogenophaga) and cluster IV (including Thiothrix, a filamentous bacteria) both displayed a decrease in quantity, whereas, cluster V had no obvious tendency. Microbial microdiversity demonstrates shifts in the dominance of functional microorganisms in ecosystems.38 To perform AOB and NOB microdiversity analyses along the RBC flowpath, sequences were clustered into species-level OTUs using QIIME, and then the representative OTU sequences (8143 OTUs) were searched for in the GreenGenes database using the Blastn algorithm with 98% sequence similarity.39 The OTUs with greater than 1.1% and 3.4% relative abundance that were classified as Nitrosomonas (3 OTUs) and Nitrospira (2 OTUs) were further examined for their relative abundance variation pattern along the flowpath. For the Nitrosomonas groups, three OTUs representing Nitrosomonas sp. Is343 and different Nitrosomonas marina, and Nitrospira (2 OTUs) comprising two groups of Nitrospira marina were analyzed. All five OTUs increased at different rates along the RBC flowpath, with the most abundant OTU of the two groups, Nitrosomonas sp. Is343 (OTU 1) and Nitrospira marina (OTU1), increasing notably along the RBC flowpath. This could be due to their tolerance to high levels of organics and ammonia, as well as to low DO (dissolved oxygen). To determine the relative increasing rate along the RBC flowpath, the relative abundances of the most abundant OTU (i.e., OTU1) were normalized to the ratio of 100 (Figure 3). The Nitrosomonas data indicated a rising trend in OTU2 (Nitrosomonas marina) in both Set A and Set B, whereas OTU3 (Nitrosomonas marina) was flattened versus OTU1. For the Nitrospira, ROTU2/ROTU1, with its ratio from 0 to 53.7%, sharply increased especially after the fifth RBC stages. The conversion of nitrogen involves multiple steps, and the increase of Nitrospira in the last RBC stage may benefit from the mutualistic cooperation between the two different functional groups. Abundance and Diversity of amoA. In this study, both the AOA and AOB amoA gene sequences were detected in all RBC stages for each train (Set A and Set B). Because AOA and AOB dominated ammonia oxidization in low ammonia environments,40 they increased sharply (both in Set A and Set B), with the lower ammonia concentrations along the RBC flowpath. This trend was consistent with a previous study on the relative abundance of RBC AOA amoA gene squences along an RBC flowpath.36 For both RBC trains, the amount of AOB amoA was greater than the AOA amoA, approximately 28 and 31 times which demonstrated a higher ratio of AOB:AOA than in activated sludge and a lower ratio than in soil.41 In addition, amoA gene relative abundances varied between Set A and Set B and were lowest in the inlet sample of Set B (RBC-B-1) and highest in the outlet sample of Set B (RBC-B-7), demonstrating that nitrifiers increased to higher ammounts in Set B than in Set A. Because Set A had been running half a year longer than Set B, biofilm biomass in Set A varied significantly with biofilm age (see Figure 4). The AOA and AOB amoA diversities across the RBC stages were analyzed by constructing a clone library. For archaeal amoA, the clone library analysis of AOA revealed that two clades were detected. Arch8 (56% of all clones, sharing 99% similarity with Crenarchaeote sp. (HQ317017)) was the most abundant, and Arch13 (18.8%, sharing 98% similarity with Thaumarchaeote sp. (KC962934)) was also found. For bacterial amoA, Nitrosomonas was the only genus detected, which was in

accordance with the community analysis from the 454 pyrosequencing data. For example, the proportion of Bact3 was 99% similar to Nitrosomonas. However, four subgroups within this genus were detected based on the amoA sequences. UniFrac Clustering of 18 Biofilm Samples. A clustering algorithm (based on Bray−Curtis distance) was applied to the 18 RBC biofilm samples to measure consistencies among various bacteria communities and to group similar biofilm samples at the genus level. As evident from Figure 5a, the

Figure 5. Cluster analysis based on the Bray−Curtis distance (a) at the genus level with a 3% cutoff. The genera were divided into a total of four groups. The four-dimensional principal coordinate plot obtained using the genus results of 18 biofilm samples is shown in (b). The samples cluster by RBC treatment stage.

bacterial communities in the 18 biofilm samples could be clustered into four groups: (1) Group I contains biofilm samples collected from RBC-A-1 to RBC-A-4 and RBC-B-1 to RBC-B-4; (2) Group II contains biofilm samples collected from RBC-A-5 to RBC-A-7; (3) Group III contains biofilm samples collected from RBC-B-5 to RBC-B-6; (4) Group V is biofilm samples from RBC-B-7. It is notable to recognize that the bacterial biofilm community shifted along the RBC flowpath. To provide additional support for this analysis, a PCoA analysis was conducted on these validation data based on relative abundances at the genus level. Figure 5b showed that the biofilm samples from the first half of both of the RBC trains (including RBC-A-1 to RBC-A-4, and RBC-B-1 to RBC-B-4) clustered together (Group I), whereas the other samples did not show that pattern (dividing into Group II, III, and IV, respectively). The biofilm samples of the RBC-A-5 - RBC-A-7, RBC-B-5 - RBC-B-6, and RBC-B-7 fall into three discrete 8049

dx.doi.org/10.1021/es5017087 | Environ. Sci. Technol. 2014, 48, 8044−8052

Environmental Science & Technology

Article

clusters. Therefore, a hypothesis can be proposed that the bacterial shift in the biofilm in the observed sequences is primarily because of the RBC treatment process of pollutant removal,42 consistent with the previous CA analysis. In summary, there are at least three hypothesized causes for the grouping of all of 18 biofilm samples into the four groups: (1) the same raw wastewater influent in both Set A and Set B led to highly similar bacteria composition in the first series, (2) a gradual bacterial community shift occurred with different pollutant removal rates, and (3) the start-up time of the biofilm influenced the biofilm formation. Correlation Analysis of Bacterial Genera in RBC. Correlations between the bacterial genealogies in the 18 RBC biofilm samples were calculated using an relative abundance matrix of genus-level OTUs with 97% similarity and were visualized in a network interface to explore significant taxon cooccurrence patterns, which could aid in deciphering the structure of complex bacterial communities and potential microbial interactions along the pollutant gradient in the RBC.41 Previously, correlation networks have been successfully applied to identify bacterial correlation and interactions in marine, soil environment, and activated sludge.43−45 In this study, we resolved the networks of interconnection among the 370 genera appearing in the 18 biofilm samples by identifying the co-occurrence patterns among them. The resulting biofilm bacterial network consisted of 36 nodes (genera) and 51 edges (connections between two genera), with an average node degree connectivity of 2.83 (representing the average number of neighbors of a node), a network diameter of 8 (the longest distance between nodes) and a modularity of 0.784 (values > 0.4 suggest that the network has a modular structure).46 The average clustering coefficient of 0.104 (high values between 0 and 1 illustrating a deviation from a randomly wired network47) is much higher than that found for an Erdös-Réyni random network (0.032) of equal size, revealing a nonrandom structure of bacterial assembly along the RBC. Overall, the RBC biofilm bacterial network comprised highly connected nodes structured within different bacterial clusters that may occupy multiple different niches in the biofilm of the RBC. A small but tightly/densely connected module of nodes was corresponded to nitrogen-cycling-related members (Figure 6), primarily including ammonia-oxidizing Nitrosomonas, nitrite-oxidizing Nitrospira, and denitrifiers such as Acinetobacter, Arcobacter, Thauera, Paracoccus, and Azoarcus. Further analysis found (P-value