Exploited Application of Pyrosequencing in Microbial Diversity of Activated Sludge System of Common Effluent Treatment Plants

Microbial Communities are actively present in the Activated Sludge System. We have applied PCRbased Pyrosequencing to investigate the bacterial communities of Activated Sludge samples from different common effluent treatment plants. A total of 259K effective sequences of 16S rRNA gene V4 region were obtained from these Activated Sludge samples. These sequences revealed huge amount of operational taxonomic units (OTUs) in Activated Sludge, that is, 1183–3567 OTUs in a sludge sample, at 3% cutoff level and sequencing depth of 16 489 sequences. Clear geographical differences among the Activated Sludge samples from effluent treatment Plant No.1 and No.2 were revealed by (1) cluster analyses based on abundances of OTUs or the genus/family/order assigned by Ribosomal Database Project (RDP) and (2) the principal coordinate analyses based on OTUs abundances, RDP taxa abundances and UniFrac of OTUs and their distances. In addition to certain unique bacterial populations in each Activated Sludge sample, some genera were dominant, and core populations shared by multiple samples, including two commonly reported genera of Zoogloea and Dechloromonas, three genera not frequently reported and three genera not well described so far. Pyrosequencing analyses of multiple Activated Sludge samples in this study also revealed the minority populations that are hard to be explored by traditional molecular methods and showed that a large proportion of sequences could not be assigned to taxonomic affiliations even at the phylum/class levels.


Introduction
The growth of the world population, the development of various industries, and the use of fertilizers and pesticides in modern agriculture have overloaded not only the water resources but also the atmosphere and the soil with pollutants [Maulin et al., 2013]. Biological wastewater treatment is a multi-billion dollar industry-the largest application of biotechnology in the world. By harnessing and concentrating microorganisms in bioreactors, the beneficial activities of naturally-occurring microorganisms are accelerated, enabling removal of oxygen-depleting organics, toxics and nutrients, and preventing discharge of pathogens. Despite the environmental and economic importance of these processes, the knowledge of the microbial communities within biological wastewater treatment is incomplete, primarily because of a lack of suitable tools for their analysis. Culture-dependent methods are biased by the selection of species which obviously do not represent the real dominance structure [ rRNA clone libraries etc., have provided new insights into microbial community structure. However these methods generally are difficult to detect most of the low-abundance organisms in a wastewater microbial community, thus provide incomplete information about the microbial diversity and community within the biological WWTPs [Xia SQ et al, 2010]. Recently developed pyrosequencing technology can generate hundreds of thousands of short sequences [Margulies M et al., 2006], and significantly improve researchers' ability to investigate the low abundance microorganisms ; Park J, 2010; Peccia J, 2010; Qian PY, 2011]. This new approach has been extended to environmental systems in recent years [Fodor, 2014]. However, to our knowledge, there are rare reports using pyrosequencing technology to compare bacterial communities between various WWTPs operated at different geographic locations. More recently, Zhang et al. [15] and Xia et al. [7] have observed that some core populations of bacteria were shared by multiple activated sludge samples from different wastewater treatment plants. However they did not examine the relationships between microbial community structures and environmental variables, thus it is not clear how wastewater characteristics, operational parameters and geographic locations influence the bacterial community structures. The low sequencing depth of the traditional PCR cloning approach when compared with the vast genetic diversity present in AS systems hindered a comprehensive characterization of the microbial community structure. In this aspect, the current community analyses typically represent a mere snapshot of the dominant members, with little information on taxa with medium to low abundances. High-throughput sequencing is a promising method, as it provides enough sequencing depth to cover the complex microbial communities [Shendure J et al, 2007]. So far, it has been applied to analyze microbial communities in marine water , soil [Roesch L, 2007;Lauber CL, 2009], human hand surface [Fierer N, 2008], human distal intestine ], etc. But few studies have been conducted using this method to investigate AS, though this approach has been used to study the microbial community in the raw sewage [McLellan S, 2010] and the residual biosolids after digestion of AS . Using computational ecology tools such as the Uni Frac bdiversity metric [Lozupone CA, 2005] and principal coordinates analysis (PCoA), researchers can compare differences in microbial communities between ecosystems and along environmental gradients. With the aid of highthroughput sequencing technology and the wellestablished b-diversity analytical tool, an attempt could be made to answer some fundamental questions related to AS microbial communities. In this study, AS samples were collected from sewage treatment plants. Pyrosequencing using the 16S rRNA gene as the biomarker was conducted to examine the bacterial diversity of these AS samples, to evaluate the similarity/difference of different samples, to compare the unique dominant bacterial populations, and to identify the core populations hared by different samples. To the best of our knowledge, this study is the first application of PCR based pyrosequencing to characterize and compare multiple AS samples.

Waste Water Treatment Plant and Sampling
Biomass samples were collected from different effluent treatment plants of India. Grab samples were collected from the aerobic basin. A bucket-type sampler was used for the suspended microorganisms, whereas a sterilized scissor was used to cut off a part of the fixed synthetic mesh, after lifting the mesh, for the attached microorganisms. Samples were immediately frozen before being transported to the laboratory.

DNA Extraction
Of these samples, the same sludge sample from effluent treatment plant was divided into two aliquots, that is, A and B. Then the two subsamples were treated as independent samples from DNA extraction to pyrosequencing, to evaluate the reproducibility of the methods applied in this study. Samples of 10ml were centrifuged at 4000 r.p.m. for 10 min at 4 0 C. Pellet (200mg) of each sample was collected for DNA extraction in duplicate with the Fast DNA SPIN Kit for Soil, which was found to be the most suitable (having the lowest contamination) for the samples in this study. The duplicate DNA extracts were then merged together for the following PCR amplification.

PCR Amplification and Pyrosequencing
The forward primer 563F (50-AYTGGGYD TAAAGNG-30) at the 50-end (E. coli positions 563-578) of the V4 region (239 nucleotides) and a cocktail of four equally mixed reverse primers, that is, R1  . A 100 ml reaction system was set up for each PCR amplification, using Mighty Amp polymerase. The amplification was conducted in an i-Cycler under the following conditions: initial denaturation at 98°C for 2 min, and 28 cycles at 98°C for 15 s, 56°C for 20 s and 68°C for 30 s, and a final extension at 68°C for 10 min. Amplicon libraries were prepared by a cocktail of three independent PCR products for each sample to minimize the impact of potential earlyround PCR errors . PCR amplicons were purified with a quick-spin Kit, and concentrations were measured using spectrometry. Amplicons from different sludge samples were then mixed to achieve equal mass concentrations in the final mixture, which was sent out for pyrosequencing at Bangalore Genei, India.

Data Analysis
Initially, the raw sequences obtained from pyrosequencing were processed using the GL FLX software for sorting by the key (i.e., sequences from the suspended and attached samples), with low quality sequences discarded, and the primer sequences trimmed. In addition, sequences less than 250-bp long were discarded. The processed sequences were subjected to taxonomical classification

Post-run Analysis
All the raw reads were treated with the Pyrosequencing Pipeline Initial Process (Cole et al., 2009) of the Ribosomal Database Project (RDP), (1) to sort those exactly matching the specific barcodes into different samples, (2) to trim off the adapters, barcodes and primers using the default parameters, and (3) to remove sequences containing ambiguous 'N' or shorter than 150 bps ). The reads selected above were defined as 'raw reads' for each AS sample. Sequences were denoised using the 'pre.cluster' command in Mothur platform to remove sequences that are likely due to pyrosequencing errors Roeselers et al., 2011). PCR chimeras were filtered out using Chimera Slayer (Haas et al., 2011). The reads flagged as chimeras was extracted out and submitted to RDP. Those being assigned to any known genus with 90% confidence were merged with the non-chimera reads, to form the collection of 'effective reads' for each AS sample. Although bacteria-specific primers were applied, very small portions of unexpected archaeal sequences might be obtained . To remove these cross-talking sequences, the effective sequences of each AS sample were submitted to the RDP Classifier (Wang et al., 2007) to identify the archaeal and bacterial sequences, and filter out the archaeal sequences using a self-written Python script. The average length of all bacterial sequences without the primers was 207 bp. After the above filtration, the minimum number of selected bacterial sequences in the 15 AS samples was 16 489. To fairly compare the different Activated Sludge samples at the same sequencing depth, normalization of the sequence number was conducted by extracting the first 16 489 sequences in each sample for all the following analyses. Taxonomic classification of the bacterial sequences of each Activated Sludge sample was carried out individually, using the RDP Classifier. A bootstrap cut off of 50% suggested by the RDP was applied to assign the sequences to different taxonomy levels. The normalized sequence set of each AS sample was also individually aligned by Infernal (Nawrocki and Eddy, 2007) using the bacteria-alignment model in Align module of the RDP. By applying the Complete Linkage Clustering, sequences in each set were assigned to phylotype clusters at two cutoff levels of 3% and 6%. The cluster analysis (CA) was conducted to group the bacterial communities of different Activated Sludge samples on the basis of (1) taxonomy results obtained using the RDP Classifier (excluding those unclassified sequences), and (2) operational taxonomic units (OTUs) generated using RDP Complete Linkage Clustering from the merged pool of sequences of all the Activated Sludge samples. A sample mapping file showing the frequency of each reference taxon in different Activated Sludge samples was generated through local BLAST.

Microbial Community Diversity
As shown in Table 1, after filtering the low quality reads using the RDP Initial Process in Pyrosequencing Pipeline (PP) and trimming the adapters, barcodes and primers, there were 20 276B28 260 effective reads for the different Activated sludge samples. After denoising, filtering out chimeras, and removing the archaeal sequences, the library size of each sample was normalized to 16 489 sequences, which were the smallest among the different Activated sludge samples, to conduct the downstream analyses for different samples at the same sequencing depth. The numbers of OTUs, Chao 1 and ACE at two cutoff levels of 3% and 6% are summarized in  Figure 2. As a common hypothesis, diversity may affect the performance of the Activated sludge process. However, an accurate estimation of OTUs in an Activated sludge sample, based on DNA sequencing, has not been conducted before. In this study, based on the 247 335 effective bacterial sequences, there could be a total of 13 951 (3% cutoff) and 8493 OTUs (6%) in the 15 samples. It should be noted that the richness values were certainly affected by sequencing noise. The 16 489 selected effective bacterial sequences in each samples were assigned to different taxa levels (from genus to phylum) using the RDP Classifier at 50% threshold. Although it has been reported that the V4 region used in this study displayed the highest number of correctly classified sequences, followed by the V3 and V6 regions , quite a large portion of effective bacterial sequences in this study could not be assigned to any taxa of different level at the 50% threshold, indicating the extent of novel sequences captured by this study. Figure 1 shows that the unclassified sequence portions in the total community increased from the domain level to the genus level, and were significantly different among the samples, especially at the family and genus levels. As shown in Figure 1 (Snaidr et al., 1997), these four groups were dominant (56-86%) in bacterial communities of the differen Activated Sludge samples in this study, followed by a few other major (average abundance 41%) phyla, including Verrucomicrobia (4.2%), Chloroflexi (3.4%), Acidobacteria (3.0%) and Planctomycetes (2.4%). A few phyla, including TM7, Thermotogae, OD1, Spirochaetes, WS3, Nitrospira and Synergistetes, were only the major (abundance 41%) phyla in one of the different samples. The abundances of other phyla were o1% in all samples. Within Proteobacteria, e-Proteobacteria only occurred at very low levels (0.01-0.51%, averaging 0.08%). Except for the four samples of effluent treatment no.5, which had the a-subdivision as the most dominant class within Proteobacteria, in all other 11 samples, the b-subdivision was the most dominant Proteobacteria, followed by a-, gand d-subdivisions. This is different from results of a study using microarray (Xia et al., 2010), which showed that a-subdivision was the most abundant of the Proteobacteria. However, the findings are similar to the results of another study about soil bacteria using pyrosequencing (Roesch et al., 2007), which demonstrated that in most soils, the b-subdivision was the most abundant one within the Proteobacteria. It is interesting that in the Activated Sludge sample from the unique saline effluent treatment plant at Ankleshwar, India the a subdivision was the dominant class within the Proteobacteria, accounting for 18.5-25.7% of total bacterial effective sequences, whereas the a-Proteobacteria only averaged 10.7% in the other 13 Activated Sludge samples. This might be explained by the dominance (up to 25-50%) of the a-Proteobacteria in the marine microbial community of the surface/subeuphotic layers (Venter et al., 2004;Qin et al., 2010). In addition, the d-Proteobacteria was the most dominant class in the other two effluent treatment plants of India, but their most dominant orders in the d-Proteobacteria were different, that is, Myxococcales (11%) and Desulfobacterales (5.8%), respectively. In addition to the four classes of the Proteobacteria, other dominant (41%) shared classes included Actinobacteria, Sphingobacteria, Clostridia, Bacilli, Planctomycetacia, Verrucomicrobiae, Anaerolineae, Verrucomicrobia Subdivision3, Flavobacteria, Caldilineae, Acidobacteria_Gp4 and Acidobacteria_Gp6. Among the 744 assigned genera, 70 (accounting for 63.7% of the classified sequences) were shared by all 15 samples. A total of 235 genera were commonly shared by more than 10 AS samples, accounting for 92.7% of all classified sequences. There were 273 rare genera that only appeared in one or two samples，accounting for only 0.9% of total sequences. a Number of AS samples, which share the genera in the second column. '•Percentage of the number of shared genera in the number of total classified genera. ' Percentage of sequences of the shared genera in the total classified sequences.

Analysis of Similarity
The similarity of the Activated Sludge sludge samples was evaluated using two independent methods: CA and PCoA. Cluster analysis. As shown in Figure 2a, CA, based on abundances of orders, revealed that bacterial communities in the 15 samples could be clustered into five groups: (1) Group I contains all Activated Sludge samples from effluent treatment plant No.1 and that from No.2; (2) Group II contains the three samples from No.3; (3) Group III is sludge from No.4, which treats saline sewage; (4) Group IV is sludge from No.5, which treats sewage containing slaughterhouse wastewater; (5) Group V is sludge from No.6, which is located inside a cave. This grouping pattern was similar at the family and the genus levels (Figure 4b, 4c), but slightly changed at the class levels (Figure 4a Figure 3 (UniFrac at 3% cut off) and Figure 5a (order) and Figure 5b (OTUs at 3% cutoff). PCoA at the family, genus and 6% cu toff levels were also conducted. Although there are slight variances among these PCoA results, the same general trend was observed, that is, the sludge samples from ETPs forming clusters different from the sludge from other ETPs. We hypothesized that there would be significant differences among the sludge from different geographical areas. As demonstrated by the CA and PCoA, the sludge samples were certainly similar to those from the same geographical area, possibly due to unique sewage compositions, as well as different plant design and operation in each area.

Shared and Distinct Orders/Families
At the order level, the 27 abundant (41% in any Activated Sludge sample) orders accounted for 45-80% of the classified sequences, as shown in Figure 3.  Figure 5 shows the genus profiles for the fifteen Activated Sludge samples. Each panel represents one of the 739 genera, which were sorted alphabetically by phyla and then genus. Examination of these genera captured the similarities and differences in the genus profiles of the 15 samples. Comparative analysis revealed a core microbiota across the 15 Activated Sludge samples. Among the 744 assigned genera, 70 (accounting for 63.7% of the classified sequences) were shared by all 15 samples. A total of 235 genera were commonly shared by more than 10 Activated Sludge samples, accounting for 92.7% of all classified sequences. There were 271 rare genera that only appeared in one or two samples, accounting for only 0.9% of total classified sequences, a very minor part of the bacterial communities in Activated Sludge. The top 10 abundant genera in each sample were selected (a total of 58 genera for all 15 samples) and compared with their abundances in other samples, as shown in Figure 4. Eight genera were abundant (41%) in at least six samples, including two genera extensively reported in Activated Sludge, that is, Zoogloea and Dechloromonas, three genera rarely reported before, that is, Caldilinea, Tricoccus and Prosthecobacter, plus three not well described genera, that is, Gp4, Gp6 and Subdivision3_genera_incertae_sedis. The species in the genus Zoogloea, such as Zoogloea ramigera, are known to form characteristic cell aggregates embedded in extracellular gelatinous matrices, often called zoogloeal matrices (Dugan et al., 1992), and are the main agent for the flocculation of Activated Sludge (Rossello-Mora et al., 1995). Caldilinea has some filamentous species and has a role in stabilizing flocs of Activated Sludge in a wide range of STPs (Yoon et al., 2010). The information about the existence and roles of other genera in AS are limited. Just like at the family level, the dominant genera showed some geographical characteristics. For example, Flavobacterium were only dominant in the three samples from Plant No.2, with abundance levels ranging from 1.83 to 7.44%, whereas its abundances were o1% in all the samples from Plant No. 3 & 4. In addition, all three samples had high levels of the photosysthetic bacteria Rhodobacter (1.43-3.72%), whereas the samples from plant no. 3 & 4 contained relatively less (0.32-0.99%). As shown in Figure 5, a few other genera more dominant included Acinetobacter, Rhodoferax, Acidovorax, Aquabacterium and Alkanindiges, whereas Gp4, Gp6, Clostridium, Hyphomicrobium, Comamonas and Pirellula were more abundant. The distribution of some abundant genera also depended on temperature. For instance, the species in Trichococcus genus have been identified as psychrotolerant mesophiles and are able to grow at low temperatures. The sludge samples also contained some other unique populations, possibly due to salinity, exposure to sunlight, industrial wastewater contribution and so on. Two of the samples from different effluent treatment plant, which had been exposed to seawater, contained unique genera such as Hoeflea, of which all three known species are related to marine environments (Palacios et al., 2006). The unique populations in different Activated Sludge samples could be due to many possible causes, including geographical isolation, concentration and types of organic substrates in the sewage, which are significantly affected by the diets of eastern and western peoples, operation mode (A/A/O or A/O vs conventional AS), dissolved oxygen concentration, temperature, salinity, season, type and percentage of industrial wastewater, etc. Although it is out of the scope of the present study, the factors shaping microbial community deserve more comprehensive and systematic studies in the future, using the methodology demonstrated here.