Sewage Reflects the Microbiomes of Human Populations

The 15 most abundant bacterial families, which on average accounted for 98% of the reads in the human stool data set, represented 26% of the sequence reads in a sewage sample (Table 1). The low representation of human fecal bacteria in sewage concurs with previous reports that 80 to 90% of bacterial sequences in sewage originate from non-human-fecal sources (24, 25). Since sewage contains many organisms of nonfecal origin, our analysis focused on sequences from bacterial families that each represented >3% of total reads in the human stool data set: Bacteroidaceae, Ruminococcaceae, Lachnospiraceae, Porphyromonadaceae, Rikenellaceae, and Prevotellaceae. Five of these represented the most abundant of the human fecal matter-associated families in sewage, with the sixth, Rikenellaceae, contributing less than Veillonellaceae (Fig. 1A). Normalization of the sewage data set to the 15 most abundant families in the human stool data set (Fig. 1A) showed comparable community compositions at the family level (linear r² = 0.72, P < 0.001) but with an underrepresentation of Bacteroidaceae (Wilcoxon rank sum test, P < 0.001) and an overrepresentation of Lachnospiraceae and Prevotellaceae (Wilcoxon rank sum test, P < 0.001).

Analysis of reads that map to the human stool sample's six most abundant families identified 351 oligotypes across all human stool and sewage samples. BLAST analyses assigned 105 oligotypes to the nonfecal data set (i.e., a set of oligotypes representing 16S rRNA genes not detected in human fecal samples; see Materials and Methods for classification details). The remaining 246 oligotypes comprised the final human fecal data set and represented on average 78% of all sequences in a human stool sample and 12% of sequences in a sewage sample (Table 1). Seven of the 246 human fecal oligotypes were detected only in the human stool samples, and nine were detected only in the sewage samples.

After filtering out the nonfecal data set sequences, oligotyping revealed an overrepresentation of Lachnospiraceae and Prevotellaceae while Bacteroidaceae remained underrepresented in sewage compared to the human stool data set (Fig. 1B). Although there existed a broad family-level composition shift between the data sets, the oligotype relative abundance profile of the pooled sewage data set resembled the pooled human stool data set (linear regression r² = 0.66, P < 0.001; see Fig. S1 in the supplemental material). Within families, oligotype relative abundance profiles correlated linearly (P < 0.001), although Lachnospiraceae had a much lower regression fit (r² = 0.42) than the other families (range from r² = 0.76 to 0.99; see Fig. S1).

Sewage exhibited the additive effect of combining microbial communities from many individuals into a single sample. Oligotypes that were common among individual stool samples were common among sewage samples (Fig. 2A) and typically more abundant (see Fig. S2 in the supplemental material). We also observed this additive effect in the alpha diversity measures for each data set. While the pooled human stool and pooled sewage data sets captured approximately equivalent oligotype richness (n = 235 for pooled stool and n = 239 for pooled sewage) and diversity (at several alpha values of the Hill diversity index [Fig. 2B]), the mean diversity for individual sewage samples exceeded values for individual stool samples (Wilcoxon rank sum test, P < 0.001 for all Hill diversity alpha values [Fig. 2B]). The high ratio for sample mean to pooled sample diversity in the sewage data set indicates that sewage had low oligotype turnover among samples relative to the low ratio (i.e., high turnover of oligotypes) among stool samples (Fig. 2B).

Sewage samples from different U.S. cities displayed very similar human fecal oligotype compositions within and among sample periods (mean Bray-Curtis similarity = 73.2% [Fig. 3]). The oligotype composition of the pooled stool data set was strikingly similar to individual sewage samples (mean Bray-Curtis similarity = 81.0%). In contrast, the mean Bray-Curtis similarity for comparisons between individual stool samples hovered at 32.6% (Fig. 3). Hill diversity measures at alpha values of ≥1, which emphasize the contribution of more abundant taxa, indicated that each sewage sample contained the majority of the most abundant oligotypes (Fig. 2B). For example, on average, two sewage samples captured 90% of the measured fecal oligotype diversity in the sewage data set, whereas the same diversity levels required data from 71 human stool samples (see Fig. S3 in the supplemental material).

No oligotype occurred in all human stool samples. The two most prevalent oligotypes were detected in 129/137 stool samples (see Table S2 in the supplemental material). In contrast, 17 oligotypes were present in all 207 sewage samples and 10 others were present in ≥205 (≥99%) of sewage samples. These 27 oligotypes also represented the most abundant amplicon sequences in the sewage samples (see Fig. S4). The 27 common and abundant oligotypes represent “core” gut microbiota among human populations in the United States. In the human stool data set, one of the core oligotypes represented the most abundant oligotype in 117 of the 137 samples. Twenty-six of the 27 core oligotypes derive from amplicon sequences that match exactly a cultivar 16S rRNA gene, and all but three oligotypes matched strains from only a single genus (see Table S2). In nearly all sewage samples (200/207), the most abundant oligotype resolved to either a Prevotellaceae, one of two Bacteroidaceae, or a Ruminococcaceae oligotype. These correspond to sequences that identically matched Prevotella sp. strain BI-42, Bacteroides dorei JCM 13471, Bacteroides vulgatus ATCC 8482, and Faecalibacterium prausnitzii A2-165, respectively.

The human fecal oligotype composition in sewage reflected increased representation of oligotypes from one family group over the others (Bacteroidaceae versus Prevotellaceae versus Lachnospiraceae plus Ruminococcaceae [see Fig. S5 in the supplemental material]; Adonis R² = 0.30, P < 0.01). Of 51 treatment plant sites that had data for all three sample collections, 21 exhibited the same enrichment patterns, a result which exceeds random expectations (binomial P ≤ 0.01). Cities with more than one treatment plant exhibited a higher level of consistency for paired sample comparisons between plants (14/19 paired samples were enriched for the same family; binomial P ≤ 0.01). Principal coordinate analyses (PCoAs) did not indicate significant (P ≤ 0.01) relationships between the geographic location of treatment plants, the plant chemical/physical measurements, or the population size served and the human fecal oligotype composition (see Table S3 for PCoA correlates). Only sample period (r² = 0.09) and city population percent obesity (r² = 0.06) explained significantly (P < 0.01) the variation in composition among U.S. cities (see Table S3 for PCoA correlates).

In contrast to the human fecal oligotypes, the nonfecal oligotype data set exhibited greater variation among cities and sample periods, with strong geographic and seasonal trends corresponding to air temperature and latitude differences. PCoAs revealed that the temperature profile of a city explained a significant (P < 0.001) proportion of the nonfecal community variation: August (r² = 0.29), January (r² = 0.36), and May (r² = 0.52). A constrained ordination for the city temperature profile parameter further illustrated that the nonfecal oligotype composition was more strongly related to a city's yearly temperature profile than the human fecal oligotype composition, and this relationship in the nonfecal community appeared as a significant divide between the northern and southern U.S. cities (Fig. 4).

Although no measured factors explained a high percentage of the human fecal community variation among cities, the percent obesity in a city's population had explanatory power: PCoA obesity percent, r² = 0.06 in August, r² = 0.08 in January, and r² = 0.11 in May (see Table S3 in the supplemental material). A random forest classification algorithm demonstrated that human fecal oligotype composition in sewage predicted whether a sample derived from a lean or obese population with 81 to 89% accuracy (Table 2). This relationship was driven partly by an increase in the relative abundance of Bacteroidaceae oligotypes in samples from the most obese city populations [see Fig. S6; Wilcoxon rank sum test for relative abundance in samples categorized as lean versus obese, W (n₁ = 51 n₂ = 54) = 806, P < 0.001], as represented by an increased ratio of two core Bacteroidaceae oligotypes versus two core Ruminococcaceae oligotypes (see Fig. S6).

Sewage influent samples were collected from 71 cities and 78 wastewater treatment plant (WWTP) sites from across the United States during August 2012 (7 August 2012 to 7 September 2012), January 2013 (9 January 2013 to 28 February 2013), and May 2013 (28 April 2013 to 4 June 2013). To obtain these samples, we shipped sampling supplies to each of the WWTPs, including a cooler, frozen cold packs, sterile 500-ml sample bottles, chain-of-custody forms, and sample instructions. WWTP operators at each site collected sewage influent according to their plant's standard collection procedures, which ranged from single-time-point grab samples to flow-weighted composite samples taken over a 24-h period (see Table S1 in the supplemental material for metadata details). Following sample collection, the operators transferred the influent to an autoclaved sample bottle and then placed the sample bottle in a refrigerator until shipment to our lab. Prior to overnight shipping, the operators sealed the sample with Parafilm and placed the sample in the provided cooler containing frozen cold packs. Upon receiving sewage influent from the WWTPs, we mixed the influent by shaking and collected the microbial communities by filtration of 10 25-ml subsamples onto 0.22-µm mixed cellulose ester filters (47-mm diameter; Millipore, Billerica, MA, USA). Each filter was stored in a 2-ml screw-cap freezer tube at −80°C for up to 3 months before extraction of DNA. Using a sterile spatula, we crushed the frozen filters in their 2-ml storage tube and added a bead-beating matrix plus buffers according to the standard protocol for the Fast DNA spin kit for soil (MP Biomedicals, Solon, OH, USA). Following bead beating for 1 min, we extracted DNA according to the manufacturer's instructions and purified sample DNA using the Mo Bio PowerClean DNA cleanup kit (Mo Bio Laboratories Inc., Carlsbad, CA, USA).

WWTP operators employed their routine procedures to collect separate samples for influent chemical/physical measurements (see Table S1 in the supplemental material for metadata details). In addition to the environmental parameters reported by the treatment plants, we aggregated data for the past 30-year average daily high and low temperatures on the collection date for each sample (data from NOAA National Climatic Data Center [http://www.ncdc.noaa.gov]), the latitude and longitude of each treatment plant, the 2010 median age for the county/counties in which each city lies (http://www.census.gov/2010census/data), and the 2009 age-adjusted percent obesity of the population (obesity considered body mass index of ≥30) for the county/counties in which each city lies (http://www.cdc.gov/obesity/data). For those wastewater treatment plants receiving water from more than one county, we calculated the average measurements of the counties involved (see Table S1 for metadata).

In total, we collected 219 samples at the WWTP sites. Fifty-seven of the 71 cities returned samples in all three sampling periods. Denver, CO; Honolulu, HI; Junction City, KS; Juneau, AK; Milwaukee, WI; Santa Barbara, CA; and Vancouver, WA had multiple sample sites within each city, and Palo Alto, CA, collected both grab and composite samples. In one case, the Fall River, MA, sample for the January period, a sample was collected outside the listed sample collection range (see Table S1 in the supplemental material), but for classification this sample was considered part of the January sampling period.

We amplified and sequenced the V4-V5 region of 16S rRNA genes (Escherichia coli positions 518 to 926, ~408 nucleotides [nt]) using primers and adaptors as described by Huse et al. (33), except for the use of a 30-cycle one-step amplification procedure with the fusion primers instead of the reported two-step amplification procedure. After size selection and quantification of the PCR products, 250 cycles on an Illumina MiSeq produced paired-end reads. To analyze the raw paired-end reads, we used the “merge-illumina-pairs” script distributed in the Illumina Utilities library (available from https://github.com/meren/illumina-utils [34]). The script removed any read-pair with more than three mismatches in the ~80-nt-long overlapping region and also required at least 66% of the nucleotides in the nonoverlapping region to have greater than a Q30 score (35). The program UCHIME (36) removed chimeras according to reference 33. After eliminating samples that did not pass quality controls at the filtration, DNA extraction, or sequencing steps, we were left with 207 sewage influent samples that were represented by 16,895,573 sequences assembled from the paired-end reads. Community sequence profiles were obtained in all three sample periods for 51 of the WWTP sites. Trimmed and quality-filtered sequences under the project name SLM_NIH2_Bv4v5 can be retrieved from the website Visualization and Analysis of Microbial Population Structures (VAMPS; https://vamps.mbl.edu/) (37).

For comparisons between sewage influent and human gut communities, we used human stool sample data from the Human Microbiome Project (5), which sequenced amplicons from the 16S rRNA gene V3-V5 region using 454 pyrosequencing technology (see reference 38 for sample collection and processing and reference 39 for sequencing-related procedures). The algorithm Global Alignment for Sequence Taxonomy (GAST) (40) assigned taxonomy to each unique stool sample V3-V5 read and each unique sewage sample V4-V5 read. We eliminated all samples that contained fewer than 4,000 sequence reads and removed one stool sample where 98% of reads resolved to the genus Propionibacterium, a nontypical genus for human gut communities (4, 7, 26). Following sample control procedures, we retained data from 137 individual stool samples (see Table S2 in the supplemental material for sample list).

The use of different primers for amplification and sequencing of the 16S rRNA region V3-V5 for the human stool data set on a 454 platform versus V4-V5 for the sewage data set on an Illumina platform could contribute to observed variation in the distribution of taxa/oligotypes between but not within data sets. To compare oligotypes from these overlapping rRNA regions, we created a combined amplicon data set for each of the six most abundant bacterial families in human stool (Bacteroidaceae, Lachnospiraceae, Ruminococcaceae, Porphyromonadaceae, Rikenellaceae, and Prevotellaceae). Following amplicon alignment with the align.seqs command in mothur (41), the nonoverlapping ends were trimmed (sequence length range after trimming, 225 to 240 nt) to create alignment files with sequences of equal length (23). A high-resolution oligotype analysis was conducted with the trimmed alignments as described previously (see reference 23 for details; oligotyping.org). Oligotyping is a supervised computational method that uses position-based Shannon entropy calculations in input alignments to identify highly variable alignment positions. The nucleotide compositions of these highly variable positions are used to parse the data into groups having identical sequences at the defined positions, and these groups are known as oligotypes (23).

To mitigate the potential biases of sampling depth on oligotyping (23), we randomly subsampled each sewage and human stool data set to a maximum depth of 80,000 and 7,000 reads per sample, respectively. The mean number of reads per sample after subsampling was 54,603 for the sewage data set and 5,996 for the human stool data set. Table S2 in the supplemental material reports the number of reads for each sample. Downstream analyses used the subsampled data sets. Oligotyping analyses minimized the impact of sequencing errors by employing a minimum substantive abundance criterion (M) and a minimum sample criterion (s). Using these criteria, oligotypes had to occur in ≥5 of 344 total samples (s) and be present at the lesser value of ≥0.01% of a family's total sequences or 500 total occurrences (M) to be included in comparative analyses. These noise-filtering steps, used previously for combined data sets (42), removed 12.5% of the human stool sequences and 2.9% of the sewage influent sequences (Table 1).

BLAST analysis against NCBI's nr database (returning the top 500 matches; July 2014) identified the highest-identity matches for each oligotype's representative sequence (unique sequence with highest count). All NCBI sequences with the highest-identity match were binned based on their NCBI database “isolation_source” record as either of human stool/gut origin or of non-stool/gut origin. We then categorized each oligotype as “human fecal” or “nonfecal” based on which bin had the highest count. For our categorization purposes and subsequent notation in the text, “nonfecal” refers to all oligotypes assigned to nonhuman gut/stool origins (i.e., these oligotypes may have highest NCBI database matches to sequences originating from fecal sources other than human). Multiple database sequence matches from the same sample were counted only once.

Oligotype categorization revealed that 105 of 351 oligotypes had best matches in NCBI's nr database to sequences of organisms or environmental samples collected from nonfecal origins (see Table S2 in the supplemental material). The remaining 246 oligotypes comprised the human fecal oligotype data set. All comparisons of sewage sequence data to the human stool samples were conducted using the human fecal oligotype data set. BLAST analyses to identify cultured isolate matches to each oligotype's representative sequence were carried out in the Ribosomal Database Project (9) sequence match program (July 2014), with “sequences ≥1200 nt” and “isolates” as exclusionary criteria.

Statistical analyses were conducted in the statistical package R (43). For all visualizations, we used the ggplot2 package (44). To calculate bacterial community similarity among samples/data sets, we used the function vegdist with the Bray-Curtis metric in the vegan package (45). Simple linear regressions were fitted with the lm function, and diversity measures were calculated with the renyiaccum function in the vegan package (45). Wilcoxon rank sum tests were used to determine whether the means between two conditions differed significantly. Heat maps were constructed using the heatmap.2 function in the gplots package (46).

We examined the relationship between the measured environmental, geographic, and demographic variables and the oligotype composition among samples with principal coordinate analysis (PCoA) using the function cmdscale in R (43). Given that temperature, latitude, and sample period are not independent factors, we ran PCoAs for each sample period in addition to analyses for all periods combined. We also created a metric that we termed the “city temperature profile” to incorporate the potential effects from climatic differences (i.e., temperature and geographic location) among cities. To calculate the city temperature profile metric, we used the previous 30-year average high and low temperatures on each of the three collection days for each city. If a city was not represented by a sample during a collection period, then the 15th day of that collection month was used in the calculation. The profile of average high and low temperatures from each of three collection days was used to construct a nonmetric multidimensional scaling (NMDS) plot (metaMDS; vegan package [45]), and each city's x axis score was used to represent the temperature profile for that city.

In order to compare community composition variation among cities between the human fecal and nonfecal oligotypes in sewage samples, we created a dummy sample matrix, where each original sewage sample was split into two samples, one containing the human fecal oligotypes and one containing the nonfecal oligotypes. The split samples were mutually exclusive for oligotypes. To account for this, each oligotype in one sample set (e.g., human fecal) was given a relative abundance of “0” in the samples in the opposing set (e.g., nonfecal). The relationship of both the human fecal and nonfecal bacterial community in each sewage sample to its city's climate was then examined in a constrained ordination of principal coordinates (CAP; capscale function, vegan package [45]) for the city temperature profile metric. The dummy relative abundance matrix allowed for the two data sets to be plotted in the same ordination and thus have comparable axis scores.

We examined the ability of the sewage bacterial community composition (human fecal oligotypes) to predict whether human populations contained a low or high percentage of obese individuals. To conduct these analyses, we created a classification model using random forests (47) implemented in randomForest v. 4.6-7 in R (48). Random forest analysis generates unbiased, “out-of-bag” error estimates for the data set without requiring the data to be split into training and test data sets. For the random forest analyses, we split the city sewage samples into three groups according to the distribution of estimated percent obese individuals for the city populations. We considered the most lean and obese human populations to be those in the first and fourth quartiles for obesity percent in our data set, which operationally resulted in a “lean” category, where ≤22.8% of individuals in a city were obese, and an “obese” category, where ≥30.4% of individuals were obese. The samples classified in the lean and obese categories were then used in the random forest model. We also compiled a second, more stringent data set using only samples that were greater than 1 standard deviation from the sewage data set's mean percent obesity. This configuration generated a lean population group consisting of samples with ≤21.5% obesity and an obese population group consisting of samples with ≥31.3% obesity.

The National Center for Biotechnology Information Sequence Read Archive has archived the raw data under BioProjects PRJNA261344 and PRJNA264400.