Molecular Diversity of a North Carolina Wastewater Treatment Plant as Revealed by Pyrosequencing

We report the results of pyrosequencing of DNA collected fromthe activated sludge basin of a wastewater treatment plant inCharlotte, NC. Using the 454-FLX technology, we generated 378,601sequences with an average read length of 250.4 bp. Running the454 assembly algorithm over our sequences yielded very poorassembly, with only 0.3% of our sequences participating in assemblyof significant contigs. Of the 117 contigs greater than 500bp long that were assembled, the most common annotations wereto transposases and hypothetical proteins. Comparing our sequencesto known microbial genomes showed nonspecific recruitment, indicatingthat previously described taxa are only distantly related tothe most abundant microbes in this treatment plant. A comparisonof proteins generated by translating our sequence set to translationsof other sequenced microbiomes shows a distinct metabolic profilefor activated sludge with high counts for genes involved inmetabolism of aromatic compounds and low counts for genes involvedin photosynthesis. Taken together, these data document the substantiallevels of microbial diversity within activated sludge and furtherestablish the great utility of pyrosequencing for investigatingdiversity in complex ecosystems.

INTRODUCTION

Top

INTRODUCTION

Although largely invisible in the urban landscape when theyare functioning well, wastewater treatment plants are integralto the municipal obligation to protect public health, aquaticecosystems, and the quality of life. At the heart of wastewatertreatment plants is a process whereby a dense microbial consortiumis employed to remove organic and nutrient contaminants. Themicrobes used to treat wastewater are a crucial tool in environmentalprotection. The current use of molecular techniques that donot require the isolation and cultivation of microorganisms(1, 33), including 16S rRNA analysis (6, 13, 20) and fluorescentin situ hybridization (8), have greatly expanded our understandingof wastewater microbial communities. Researchers have identifiedmany bacteria of importance to wastewater treatment, includingthe bacteria involved in biological phosphorus removal (5, 16,29), nitrifiers (8, 19, 25), denitrifiers (3, 12, 17), and methanogens(18, 36). Molecular techniques have also improved our understandingof fundamental processes such as nitrification and denitrification,as well as plant upsets, such as foaming (9, 24), which candecrease treatment efficiency.

In this paper, we apply recently developed pyrosequencing technologyto probe the molecular diversity of the aerobic basin of a wastewatertreatment plant in Charlotte, NC. In line with other studiesof complex microbial communities (28, 32), we observed astoundinglevels of diversity. We found that substantial regions of thegenomes of the most prevalent microbes in the wastewater treatmentplant are poorly described by existing sequence databases. Ourresults demonstrate that despite recent technological advancesthat allow identification of microorganisms, the microbial populationof wastewater treatment plants remains undersampled and inadequatelycharacterized. Our results are a first step toward more completemolecular characterization of this important microbial community.

MATERIALS AND METHODS

Top

MATERIALS AND METHODS

The Mallard Creek Water Reclamation Facility is located in Charlotte,NC. The plant has an average daily inflow of 7.5 million gallons,and the wastewater is mostly domestic, with additional inputfrom the University of North Carolina at Charlotte, UniversityCity Carolinas Medical Center hospital, and several industrialusers. A schematic diagram of the flow through the plant isshown in Fig. S1 in the supplemental material. Influent rawwastewater is screened and sent through grit removal beforeit is routed to day tank equalization basins that distributethe flow among three primary clarifiers. Primary effluent entersanoxic basins, where it is joined by recycle flow from the aerationbasins. Effluent from the anoxic basins enters aeration basins(solids retention time, $~$ 8 days) and then flows to secondaryclarifiers. Clarified effluent is routed to denitrificationfilters and then to UV disinfection before discharge into MallardCreek.

The plant National Pollutant Discharge Elimination System permitrequires the plant to meet a monthly 5-day test for carbonaceousbiochemical oxygen demand of 4.2 mg/liter in the summer and8.3 mg/liter in the winter months. Ammonia nitrogen (NH₃-N)levels must be below 1 and 2 mg/liter in the summer and winter,respectively. There are no other nitrogen or phosphorus limits.The total suspended solids are limited to a maximum of 30 mg/liter,and the pH must be between 6 and 9. Fecal coliform counts mustbe less than 200 CFU per 100-ml sample. These limits are routinelymet by the plant unless there are extreme weather events orplant upsets. Wastewater entering the secondary treatment systemwas monitored over a 6-month period for filtered flocculatedchemical oxygen demand, a good estimator of readily biodegradablesoluble organics, and the values ranged from 40 to 75 mg/liter.The ammonia nitrogen concentrations in this same flow rangedfrom 12 to 24 mg/liter, with the concentration varying in partdue to return flow from digested sludge dewatering.

On the morning of 20 March 2007 we collected a 50-ml samplefrom the aeration basin using a plastic dipper. At the timeof sample collection, the temperature in the aeration basinwas 18.5°C and the pH was 6.5. The sample was decanted toremove as much foam as possible before the liquid was transferredto a sterile tube. DNA was extracted from the sample using aMo Bio UltraClean Water DNA kit. The sample tube was invertedseveral times to maximize homogeneity, and a 10-ml aliquot wasremoved and pipetted onto the provided filter (0.22 µm).The filtrate was discarded, and DNA was extracted from the membraneusing the manufacturer's protocol. The final DNA extract wasanalyzed to determine its purity and concentration using a NanoDropND-1000 spectrophotometer. Approximately 100 µl of extractedDNA was concentrated in a vacuum centrifuge and resuspendedin about 12 µl of molecular-grade biology water. The finalsample concentration was 479 ng/µl as determined by aNanoDrop spectrophotometer. Preliminary analysis of the DNAusing denaturing gradient gel electrophoresis indicated thatthere was substantial diversity in the observed bands, confirmingthat our DNA extraction was successful (data not shown). Thesample was submitted to 454 Life Sciences for pyrosequencingby the 454-FLX platform. The methodology underlying pyrosequencinghas been documented elsewhere (22).

The bioinformatics analyses used in this study are shown inFig. S2 of and described in the supplemental material.

Nucleotide sequence and quality score accession number.
Sequences and quality scores from our pyrosequencing run havebeen deposited in the NCBI short-read archive under accessionnumber SRA001012.

RESULTS AND DISCUSSION

Top

RESULTS AND DISCUSSION

Our sequence set largely fails to assemble, although contigs that were generated from the assembly include many transposons as well as many genes encoding hypothetical proteins.
Our pyrosequencing run yielded 378,601 sequences with a readlength of 250.4 ± 29.1 (mean ± standard deviation).The distribution of sequence lengths was approximately normal,with a small left tail indicating some short reads (see Fig.S3 in the supplemental material). We attempted to assemble sequencesin this data set using version 1.1.02 of the GS De Novo Assemblerof the Genome Sequencer FLX data analysis suite with the defaultparameters applied. This assembly algorithm attempts to combineindividual sequence reads into longer "contigs." Given thatmetagenomic data sets of complex ecosystems have been extremelyresistant to assembly (23, 28), we expected to see very littleassembly in our data set. The 454 sequence assembler definesa "large" contig as a contig that consists of at least 500 bp.Because our average sequence length was $~$ 250 bp, this thresholdcould be achieved with overlap of a modest number of our sequences.Despite this, only 1,154 (or approximately 0.3%) of our readswere recruited into 117 contigs greater than 500 bp long (thesequences of these contigs are available in File S1 in the supplementalmaterial). To assign possible functions to these contigs, weused the GenMark algorithm (4) to predict genes on our contigsand then performed a Blastp search of these predicted proteinsagainst the pfam database. This method produces more assignmentsthan other approaches, including those based on profile searches(see the supplemental bioinformatics methods for details). Withan E-value cutoff of 0.01, this approach found matches for 75%(88/117) of our large contigs (see File S2 in the supplementalmaterial). Of these matches, 22% (20/88) were to hypotheticalproteins and 21% (19/88) were to transposases. The prevalenceof transposase gene sequences in our assembled contigs suggeststhat transposons are much more strongly conserved across metagenomesthan other genomic regions. The prevalence of contigs for hypotheticalproteins shows that the function of many of the highly conservedregions of our metagenome is poorly understood.

This failure of the 454 assembly algorithm to assemble 99.7%of our sequence reads emphasizes the great diversity of themicrobial community within the treatment plant. Because previousstudies found a similar failure of assembly algorithms for metagenomiccommunities characterized by Sanger sequencing (28), as wellas for simulated data sets created by sampling Sanger sequencingreads (23), we would not expect a significantly improved degreeof assembly even if our sequence reads were longer.

The majority of taxa in the wastewater treatment plant cannot be classified at the genus level.
In order to discover the 16S rRNA genes within our data set,we downloaded the 16S rRNA gene FASTA DNA sequences from version9.52 of the Ribosomal Database Project (RDP) (7) and used thesesequences to create a BLAST database. Using the Blastn algorithm,we asked which of our 378,601 query sequences could be foundin this RDP database with an E-value of $≤$ 0.01 (see supplementalbioinformatics methods for details). The resulting 648 sequences(available as File S3 in the supplemental material) were runthrough the RDP classification algorithm (34). The RDP classifieralgorithm uses Bayesian statistics to assign taxa to 16S rRNAgene sequences. The output of this algorithm includes a confidencescore, which ranges from 0 to 100, that indicates the degreeof confidence that can be assigned to the classification basedon the results of 100 bootstrap trials (see reference 34 formore details). The recommended threshold for assignment of ataxon by the RDP algorithm is a confidence score of $≥$ 80. Becausesequence reads as short as 90 bp have been shown to be longenough to accurately characterize taxa (15, 21), we anticipatedthat our results would not be substantially different even ifwe had a read length of more than 250 bp.

The classifications of the 148 16S rRNA sequences that couldbe assigned to a phylum with a confidence score of $≥$ 80 are shownin File S4 in the supplemental material and are summarized inFig. 1. In another paper (T. J. Hamp, W. J. Jones, and A. A.Fodor, submitted for publication), we show that these classificationsof 16S rRNA sequences derived from the whole-genome wastewatersequence set are well correlated with results from PCR experimentstargeting the 16S rRNA gene. At the phylum level, the observedtaxa are dominated by the Proteobacteria, with $~$ 70% of the classifiabletaxa belonging to this category (Fig. 1, top panel). Movingfrom phylum to genus, fewer of the sequences can be classifiedwith an RDP confidence score of at least 80%. At the genus level,nearly 60% of the sequences cannot be classified at an RDP thresholdof 80, and, among the taxa that can be classified, there isno dominant taxon (Fig. 1). These data demonstrate the extraordinarymicrobial diversity of activated sludge and are consistent withreports for other complex environments (14, 28, 30). We notethat the inability of the RDP algorithm to classify these sequencesto taxa with high confidence is not primarily because our 16SrRNA sequences have never been observed previously. Figure 2shows that many of the sequences with RDP scores of <80%(to the left of the vertical lines) have very high levels ofidentity with previously described sequences. These resultsdemonstrate that for wastewater treatment plants, as is thecase for other complex ecosystems, the accumulation of 16S rRNAsequences in public databases is vastly outpacing our abilityto classify these sequences and that this problem becomes morepronounced as one moves from the phylum level to the genus level.Presumably, future annotation efforts will rectify this problem.

View larger version (46K):
[in this window]
[in a new window]

FIG. 1. Pie charts showing taxonomic assignments for 148 16S rRNA sequences in our data set that could be classified to the phylum level with RDP confidence scores of $≥$ 80. At the phylum level, the Simpsons diversity index is 0.48.

View larger version (24K):
[in this window]
[in a new window]

FIG. 2. Results obtained with the RDP classification algorithm for 148 16S rRNA sequences that can be assigned at the phylum level with a confidence score of $≥$ 80. The x axis of each graph shows the confidence in assignments as reported by the RDP classification algorithm. The y axis of each graph shows the level of identity (expressed as a percentage) between our query sequence and the best Blastn hit in the RDP database (version 9.52). The horizontal and vertical lines indicate 95% sequence identity and an RDP confidence score of 80, respectively.

16S rRNA gene sequences from freshwater, soil, and other wastewater studies dominate our sequence set.
For each of the 648 sequences in our pyrosequencing data setthat matched sequences in the RDP 16S rRNA database (version9.52) at an E-value cutoff of $≤$ 0.01, we manually annotated wherethe corresponding RDP sequence was discovered. This was doneby manual inspection of the GenBank records for these 648 sequences.The results of this annotation are shown in File S5 in the supplementalmaterial and in Fig. 3. The x axis of Fig. 3 indicates our classification,while the y axis indicates the E-value with which the top hitfrom each of our query sequences matched the RDP database sequence.We found that while a large number of environments had at leastone hit, if we restricted ourselves to environments with multiplehits at high stringency (i.e., low E-value), only three environmentsare well represented: freshwater, soil, and other wastewaterstudies (Fig. 3). While, of course, the low number of sequencesfor some of the other environments may simply reflect the lownumber of sequences from those environments in the RDP 16S rRNAdatabase, there is a strikingly small number of sequences withhigh scores that are related to two 16S rRNA populations thatare well represented in the database: marine and human. Therelatively small number of human-derived 16S rRNA sequencesobserved is particularly interesting given the vast number ofhuman microbes deposited in the wastewater treatment plant eachday. These results show that the environment within the wastewatertreatment plant exhibits strong selection pressure against themicrobes that are present in human feces.

View larger version (12K):
[in this window]
[in a new window]

FIG. 3. Locations (as determined by manual annotation) and E-values of sequences from the 648-member pyrosequencing data set that matched the RDP 16S rRNA database at an E-value cutoff of 0.01.

Sequenced bacterial genomes are not well represented in the wastewater metagenome.
When Blastn is used to compare our sequences to the nt database,only 34% of the sequences (73,274/378,601) match the nt databaseeven at a relaxed threshold E-value cutoff of 0.01 (data notshown; see the supplemental bioinformatics methods). For thevast majority (over 98%) of the sequences that do match thent database at this threshold, the best hit is to bacterialtaxa (data not shown; see the supplemental bioinformatics methodsfor details). Since wastewater treatment plants are known toharbor many eukaryotes (27), this result likely reflects ourDNA isolation strategy, which was designed to capture prokaryoticDNA, rather than the "true" ratio of prokaryotes to eukaryotesin the treatment plant.

As of November 2008, there were 772 complete bacterial genomesin the NCBI database (ftp://ftp.ncbi.nih.gov/genomes/Bacteria/all.fna.tar.gz).In order to explore how well these known genomes are representedin the treatment plant, we used Blastn to compare our wastewatersequences to the 1,442 assembled genome and plasmid sequencesfrom the 772 sequenced bacteria. In order to eliminate spurioushits, we required that any hit matched at least 75 nucleotidesin our query sequence (see the supplemental bioinformatics methodsfor more details). Because our average sequence length was 250.4bp, this is not an overly conservative criterion. Using thiscriterion, only 20% (73,274/378,601) of our sequences matchedany of the known bacterial genomes. This result again reflectsthe great diversity of organisms in the wastewater treatmentplant and emphasizes a key challenge for genomics; despite theconsiderable effort that has been expended in microbial genomeprojects, the great majority of our sequence reads are not foundin known genomes.

For the sequences that do match known genomes, we can determinehow closely the sequenced genomes of cultivated organisms matchthe genomes present in our wastewater metagenome. We calculatedfor each of the 1,442 assembled sequences from the 772 finishedgenome projects the number of nucleotides in that genome thathave a Blastn match that aligns with at least one of our wastewatersequences. Dividing this number by the total length of eachassembled sequence yielded the "fraction of the genome covered."Figure 4 shows that even for the bacterium with the most well-representedassembled genome, nitroaromatic compound degrader Acidovoraxsp. strain JS42 (accession number NC_008782), only 25% of itsgenome sequence matched our wastewater metagenome. Table 1 showsthat the fraction of the genome covered is similarly poor forthe 10 genomes that recruited the most reads from our wastewatermetagenome.

View larger version (22K):
[in this window]
[in a new window]

FIG. 4. Fraction covered as a function of the size of each assembled sequence for each of the 1,442 assembled plasmids and chromosomes in the NCBI datadase. The fraction covered is defined as the number of nucleotides in the assembled sequence that match at least one of our wastewater sequences divided by the total number of nucleotides in the assembled sequence.

View this table:
[in this window]
[in a new window]

TABLE 1. Top 10 assembled microbial genomes as sorted by the number of hits recruited from our wastewater metagenome^a

Figure 5 shows a recruitment graph for our wastewater treatmentplant for the assembled Acidovorax sp. strain JS42 chromosome,which recruited the most sequences from our wastewater genome(Table 1; see File S6 in the supplemental material). The x axisshows the positions where sequence reads are mapped with Blastnagainst the Acidovorax genome. The y axis shows percent identityof a read compared to the matching subsection of the Acidovoraxgenome. Figure 5 shows sequences from two different sources:our 20 March aeration basin pyrosequencing run and the environmentalsequence database from NCBI downloaded in June 2007, which atthat time was largely dominated by sequences from the J. CraigVenter Institute's Global Ocean Sampling (GOS) expedition (28).We included the environmental sequence database because we wantedto assess how specific our wastewater treatment plant resultswere relative to other metagenomic sequencing databases.

View larger version (19K):
[in this window]
[in a new window]

FIG. 5. Nonspecific recruitment against the Acidovorax sp. strain JS42 genome. BLAST hits with alignment lengths less than 75 nucleotides (for the 20 March run) or 250 nucleotides (for the environmental sequence database) were removed. Protein annotations are derived from the full NCBI core nucleotide report for the Acidovorax sp. strain JS42 genome (http://www.ncbi.nlm.nih.gov/entrez/viewer.fcgi?db=nuccore&id=121592436).

The pattern shown in Fig. 5 is typical of nonspecific recruitment.For regions of the genome with conserved genes, both sourcesof sequence matched the genome, but the levels of identity wereusually below 90%. For regions of the genome that are poorlyconserved, such as the region encoding the putative transmembraneprotein (Fig. 5), very few sequences from either source mappedto the genome. We observed similar patterns of nonspecific recruitmentover a number of the genomes that recruited large numbers ofsequence reads in our 20 March aeration basin data set (datanot shown). In the Venter Institute's GOS survey, a similarpattern of nonspecific recruitment was observed for nearly everyknown microbial genome despite the presence of over 7,600,000sequences in the data set (28). This result is one of the principalreasons that the GOS study concluded that microbial diversityin the oceans is profound (28). Our results show that the genomesof the most abundant microbes in the wastewater treatment plantare largely uncharacterized. Moreover, the pattern of nonspecificrecruitment shown in Fig. 5 suggests that even additional whole-genomeshotgun sequencing would not improve the match between knowngenomes and the sequences observed in our metagenome.

The great diversity of our wastewater metagenome caused veryfew contigs to be assembled. Of the sequences that were joinedas contigs, a substantial fraction involved transposases. Wemight expect, therefore, a different pattern of recruitmentaround transposons. Figure 6 shows a region of the Acidovoraxgenome around a transposase gene sequence with a stark exceptionto the pattern of nonspecific recruitment. A large number ofsequences from our metagenome recruited to this region witha nearly perfect match. Interestingly, a number of marine sequencesfrom the GOS (28) also matched the region around this transposasegene, suggesting that, unlike most genomic regions, parts ofthis transposon are conserved across a wide environmental space.

View larger version (17K):
[in this window]
[in a new window]

FIG. 6. Region involving a transposase from the JS42 genome that shows an exception to the pattern of nonspecific recruitment. For visualization, a small amount of random noise was added to the y axis (as otherwise most of the hits to the transposase region would be superimposed). The red sequences matching the transposase region are from the GOS (28).

One genome of particular interest that has not yet been depositedas an assembled genome in the NCBI database is the genome of"Candidatus Accumulibacter phosphites," a taxon that dominatestwo recently sequenced lab-scale enhanced biological phosphorusremoval sludges (11). Although the fully assembled genome wasnot available, we did perform a BLAST comparison of our wastewatermetagenome to the largest assembled contigs from the "CandidatusAccumulibacter" project, and we observed nonspecific recruitment(data not shown). This suggests that the "Candidatus Accumulibacterphosphatis" taxon is not a dominant member of our North Carolinawastewater metagenome.

When mapped to protein space, the wastewater metagenome displays a distinct metabolic profile.
By translating our nucleotide sequences in all six frames andmapping the translated sequences to known proteins, we can generatea distinct metabolic profile for our wastewater sequences. Thisapproach, asking which genes a microbial community is capableof producing, has been successfully used to analyze the metabolicsignatures of a number of metagenomic sequence sets (10, 31).To perform this analysis, we submitted our pyrosequencing dataset for annotation on the SEED platform (2, 26). Within SEED,metabolic pathways are classified in a hierarchical structurein which all of the genes required for a specific task are arrangedinto subsystems. At the highest level of organization, the subsystemsinclude both catabolic and anabolic functions (for example,DNA metabolism), and at the lowest levels the subsystems arespecific pathways (for example, the synthesis pathway for thymidine).Using the Blastx algorithm and an E-value cutoff of 0.001, theSEED database was able to assign $~$ 60% of our sequences. The resultof assigning these sequences to functional categories is shownin Fig. 7. For comparison, Fig. 7 shows the mapping to functionalcategories from a recently published survey of 1,040,665 sequencesfrom 45 microbial metagenomes collected from nine distinct biomes(10). We note that compared to the "average" profile of thesenine biomes, the wastewater treatment plant has a distinct metabolicsignature. For example, compared to other biomes, the wastewatertreatment plant contains almost no genes coding for proteinsinvolved in photosynthesis. We would expect this as the primaryenergy source for the microbes at this treatment plant is theorganic material being processed by the plant. In addition,genes involved in the degradation of aromatic compounds areexpressed at a much higher rate in the wastewater treatmentplant than in other metagenomic systems. Again, we might expectthis given the nature of household and industrial wastes presentin sewage. Finally, we note that the Mallard Creek WastewaterTreatment Plant has no additional biological nutrient removalfacilities to treat phosphorus. Consistent with this, the percentageof sequences assigned to genes involved in phosphorus metabolismappears to be lower than that for genes involved in nitrogenmetabolism in the activated sludge (Fig. 7).

View larger version (39K):
[in this window]
[in a new window]

FIG. 7. Functional categories provided for our data set by the SEED server (http://www.theseed.org). The data for microbial genomes are averages for sequences gathered from multiple biomes (10).

Summary.
We are at the beginning of a sequencing revolution. The 91 millionbase pairs in the sequence data described in this paper weregenerated from a single sequence run with a 454-FLX instrumentgenerating over 6,000 bp of sequence per dollar. This is anapproximately 10-fold-lower cost per base pair than Sanger sequencing,and moreover, this procedure eliminates the costly and time-intensivestep of creating a bacterial clone library. As new sequencingtechnologies continue to be developed, we can expect both thecost and the experimental effort associated with metagenomicsequencing projects to drop exponentially.

Perhaps the most surprising result of our study is the pronouncedconservation of transposases across widely different environments.While there is generally poor agreement between sequences fromthe GOS and known genomes (28) and between our wastewater genomesand known genomes (Fig. 4 and 5), there are a few regions ofconservation involving transposons (Fig. 6) where there is apronounced match between the metagenomes and the sequenced genomes.A substantial fraction of the contigs that could be assembledfrom our data set involved strongly conserved transposases.It is an open question why transposons have escaped the pronouncedsequence mutability that mark nearly all of the rest of bacterialgenomes.

Like the results of other metagenomic projects (28, 32, 35),our results point to the extraordinary diversity of microbialcommunities. Patterns of nonspecific recruitment to known genomessuggest that the structures of the genomes of the most abundantorganisms in the wastewater treatment plant are unknown (Fig.4 and 5). Despite the great diversity of microbes in the treatmentplant, analysis at the protein level is surprisingly tractable,with the sequences from the treatment plant displaying a distinctmetabolic profile consistent with what we would expect basedon the plant's function (Fig. 7). This suggests that despitethe great complexity of microbial communities, next-generationsequencing technology will be a useful tool for monitoring changesin microbial processes across time and space. As treatment requirementsbecome more stringent and monitoring expands to address a broadeninggroup of compounds of concern, probe-free sequencing will increasethe rate at which key microbial groups can be identified andselected for to optimize contaminant removal.

FOOTNOTES

* Corresponding author. Mailing address: Bioinformatics Research Center, Cameron 212, UNC Charlotte, 9201 University City Blvd., Charlotte, NC 28223. Phone: (704) 687-8214. Fax: (704) 687-6610. E-mail: anthony.fodor{at}gmail.com

Published ahead of print on 29 December 2008.

Supplemental material for this article may be found at http://aem.asm.org/.

REFERENCES

Top