Metagenomics of urban sewage identifies an extensively shared antibiotic resistome in China

1Antibiotic resistance is one of the most serious global threats to human health, challenging the treatment of life-threatening infections. The widespread use of antibiotics in humans and animals is the main selective driving force of the emergence and dissemination of antibiotic resistance, and thus the cure is also the cause. Antibiotic resistant pathogens now occur at high frequencies in clinical contexts, and are increasingly being found in environmental settings, such as water bodies, soils, and animal feces. In particular, the frequent presence of multi-antibiotic resistant “superbugs” in human feces predicts a return to the pre-antibiotic era, where a growing number of infections can no longer be treated using the current arsenal of drugs.

2Consequently, the World Health Organization has endorsed a Global Action Plan on Antimicrobial Resistance, which calls upon all nations to adopt mitigation strategies within 2 years. However, there is still more to understand about the ecology and evolution of antibiotic resistance. In particular, not enough is known about the properties of the microbial resistome in ecosystems dominated by humans, and how to monitor such environments in order to evaluate their potential for promoting the evolution of antibiotic resistance.

3Municipal wastewater treatment plants (WWTPs) receive and digest millions of tons of domestic sewage. Adults harbor significant quantities of resistance genes in their gut microbiome, and consequently WWTPs, especially influents are likely to be a critical hub for the evolution and spread of anthropogenically derived resistance genes into natural environments. In China, more than 3700 municipal WWTPs have been constructed to treat urban sewage, with a combined capacity of 157 billion liters per day. In each of these facilities, sewage from tens to hundreds of thousands of individuals creates an enormous biological reactor where bacteria, and resistance genes are exposed to significant concentrations of selective agents such as antimicrobial agents, disinfectants and heavy metals. The selective pressures exerted by these agents, together with the presence of dense bacterial populations facilitates selection of antibiotic resistance and the generation of additional resistant bacteria via horizontal gene transfer (HGT). This makes sewage a vast repository of bacteria that carry and exchange resistance genes. In this respect, resistance genes detected in sewage might represent the resistance burden of their urban populations. Resistance profiles in sewage could then reflect the structure and diversity of resistant bacteria in the gastrointestinal tracts of urban residents within the WWTP catchment. This may be especially true when the WWTP mainly treats domestic wastewater without significant contributions from agricultural and industrial sources. A nation-wide survey of resistance elements in sewage (untreated influent) could then provide a rapid and efficient method for assessing the burden of antibiotic resistance from urban populations.

4Urban sewage compositions are subject to strong temporal and environmental variation in conditions. However, if and how the composition of microbial community and antibiotic resistomes change with seasons and regions in urban sewage have not been extensively investigated. To address this need, 116 urban sewage samples were collected from 32 WWTPs in 17 major Chinese cities during summer and winter. Sampling sites were specifically chosen to reflect diverse climatic conditions, economic development levels and urban geography. By combining metagenomic analyses and Illumina sequencing of 16S rRNA genes, the seasonal and geographical variations of antibiotic resistome and corresponding microbial community structure were characterized.

5A total of 116 sewage samples were collected from 32 WWTPs influents in 17 major Chinese cities during August 2014 (summer, n = 59) and February 2014 (winter, n = 57). All the untreated influent samples from each WWTP were taken within two consecutive days without recent rainfall to exclude the effect of the weather. Detailed information on these samples is summarized in Table 1 and Additional file 1: Table S1. All sewage samples from WWTPs were collected in 400-mL sterilized containers and were mixed with 100% ethanol at a volume ratio of 1:1 for biomass fixation. The fixed samples were kept on ice and were immediately delivered to laboratory. The microbial cells from 400 mL of fixed sample were pelleted by centrifuging at 9500 g for 20 min at 4 °C. All pellets were stored at −20 °C before DNA extraction.

6Genomic DNA was extracted from the collected pellets using the FastDNA® Spin kit for Soil (MP Biomedicals, France) following the manufacturer’s instructions. Total DNA was eluted in 100 μL of sterile water and kept at −20 °C until use. DNA concentrations and purity were measured using a NanoDrop spectrophotometer (ND-1000, Nanodrop, USA).

7The hypervariable V4-V5 region of the 16S rRNA gene was amplified using the primer pair (515 F: 5′-GTGCCAGCMGCCGCGG-3′ and 907R: 5′-CCGTCAATTCMTTTRAGTTT-3′ with sample-identifying six-nucleotide barcodes). The 4 × 50 μL reaction system was set up for each PCR amplification under the following program: initial denaturation at 95 °C for 5 min, and 30 cycles at 95 °C for 30 s, 58 °C for 30 s, and 72 °C for 30 s and a final extension at 72 °C for 10 min. The resulting amplicons were purified, quantified, pooled, and sequenced on an Illumina MiSeq PE300 platform (Novogene, Beijing, China). For metagenome sequencing, approximately 3 μg of sewage DNA was used for shotgun library construction with an insert size of 300 bp, followed by Illumina paired-end sequencing on the HiSeq 2500 platform (Novogene, Beijing, China).

8All the raw reads were processed using QIIME (1) to sort and assign by exactly matching the unique barcode into each sample, (2) to remove primers and the sequences with ambiguous bases, primer mismatches, and homopolymers in excess of six bases, or error in barcodes, and (3) to filter low-quality reads with >20 low-quality bases. Chimeric and noisy sequences were also filtered out. After processing, sequences were clustered into operational taxonomic units (OTUs) using Uclust clustering, which groups sequences at a minimum pair-wise identity of 97%. Mitochondrion or chloroplast sequences and singleton OTUs were discarded from the final OTU table. For each resulting OTU, the most abundant read was selected as a representative sequence. The taxonomic classification of each representative sequence was conducted using a Ribosomal Database Project (RDP) Classifier at an 80% confidence threshold (Version 2.2). Alignment of the OTU representative sequences was conducted using a PyNAST aligner, and a phylogenetic tree was built using a FastTree algorithm for downstream analysis.

9Rarefaction was performed to discern Phylogenetic Diversity, Chao1 diversity, Shannon index, and observed species metrics at each sampling depth. To remove the bias caused by different sequencing depth, the OTU table was rarefied and an even sampling depth was set by randomly subsampling the same number of sequences from each sample. Beta-diversity was estimated by computing weighed/unweighed UniFrac and Bray-Curtis distances between every pair of community samples using QIIME.

10Of the original 116 sewage samples, 24 samples were excluded from metagenomic sequencing due to the low quantities of DNA or poor sequence data (Additional file 1: Table S1). Thus, only 92 samples were further used for metagenomics analysis. Metagenomic sequencing of sewage DNA samples generated 203 Gb pairs of high-quality data with an average of 2.2 Gb for each sample. Data filtration was conducted to remove raw reads with low-quality following the methods used in a previous study. Subsequently, metagenomic sequences were analyzed by BLASTx against the Structured Non-redundant Clean Antibiotic Resistance Genes Database (SNC-ARDB) with E value ≤1 × 10−5. A read was annotated to be a resistance gene if its BLAST hit for the alignment against SNC-ARDB had ≥90% amino acid read identity for ≥25 amino acids. In the present study, a package of customized scripts was developed to automatically classify the BLAST hits into different types and subtypes of resistance genes. The detailed procedure for sorting sequences using a customized Python script was reported previously.

11SNC-ARDB contains a number of genes for efflux proteins that do not necessarily confer resistance phenotypes. These proteins do, however, function in the efflux of antibiotics and have previously been classified as resistance genes and in the Comprehensive Antibiotic Resistance Database (CARD). Therefore, efflux pump-related genes were retained in the SNC-ARDB to evaluate antibiotic resistance potential.

12The ‘abundance’ of the resistance type or subtype was calculated as previously reported by Li et al. (2015). Thus, the abundance of resistance genes based on metagenomic analysis was compared with those derived from qPCR in the previous studies. The abundance of resistance genes was transformed to ‘concentration’ (copies per liter) by normalization to the absolute copy number of 16S rRNA gene. The average copy number of 16S rRNA genes per bacterium is currently estimated at 4.1 based on the Ribosomal RNA Operon Copy Number Database (rrnDB version 4.4.4). The numbers of bacterial cells were calculated by dividing the copy number of 16S rRNA gene by 4.1, and the ‘relative abundance’ of resistance genes (copies per bacterial cell) was estimated by dividing the ARGs concentration in each sample by its corresponding number of bacterial cells. Additionally, the copy number of resistance genes discharged by each person per day in urban areas is defined as ‘ARG load’, which can be calculated by the formula: ARG load (copies/capita/day) = (the average concentration of sewage ARGs) × (the volume of municipal sewage discharge)/urban population. ‘ARG burden’ (copies/day, the total ARG load in urban areas) is calculated by multiplying the medium value of ARG load by total urban population.

13Real-time qPCR assay of total bacteria was performed using a SYBR® Green approach on a Roche 480 (Roche Inc., USA). The absolute copy numbers of 16S rRNA gene were quantified using primers 515 F and 907R. The qPCR system (20 μL) amplification was conducted as reported previously. The size of amplified fragments was about 410 bp. For the preparation of 16S rRNA gene standards, 16S rRNA gene was amplified from extracted DNA and then was cloned into the pMD 19-T vector (TaKaRa, Japan). Plasmids containing the target gene were used as standards for the calibration curve. All qPCR assays were conducted in triplicate with negative and positive controls.

14In this study, human gut microbiota was defined as the bacterial genera detected in human gut or human intestinal tract. To track the human gut microbiome fingerprint in the sewage, the human gut microbiome database including 382 bacterial genera (Additional file 2: Dataset 1) was retrieved from Human Microbiome Project (HMP) and the Metagenomics of the Human Intestinal Tract (MetaHit) project. The microbial catalogue reference set (16S rRNA gene) from these two human metagenomic projects covers almost all genera of human gut bacteria and is a useful resource for further analyses of human gut microbiome.

15Averages and standard deviations were determined using Excel 2010 (Microsoft Office 2010, Microsoft, USA). One-way analysis of variation (ANOVA), paired-sample t tests and correlation tests were performed using SPSS V20.0 (IBM, USA). All statistical tests were considered significant at P < 0.05. Diversity index, non-metric multidimensional scaling (NMDS) and significance test (Adonis test, procrustes analysis, and mantel test) were performed in R 3.1.0 with vegan 2.2.0. Post-hoc plot was generated using STAMP V2.1.3. To investigate co-occurrence patterns of microbial community and resistome, correlation matrixes were constructed by calculating each pairwise Spearman’s rank correlations. The P value was adjusted with a multiple testing correction using FDR method to reduce the false-positive results. A correlation between any two items was considered statistically robust if the Spearman’s correlation coefficient (ρ) was > 0.7 and the P value was < 0.01. The resulting correlation matrixes were translated into an association network using Gephi 0.9.1. An informatics mathematical approach based on geographical information systems, ArcGIS, was applied to map the resistance load and the resistance burden at varying spatial scales. Spatial autocorrelation analysis was conducted to evaluate spatial dependency of ARG burdens between provinces using ArcGIS.

16Twenty resistance gene types consisting of 381 subtypes were detected, with 373 subtypes in summer samples and 346 subtypes in winter samples, respectively (Fig. 1a). The three most dominant resistance gene types, conferring aminoglycoside, tetracycline, and beta-lactam resistance, accounted for 54.1% of the total ARG abundance (Additional file 3: Figure S1a). For the resistance gene subtypes, genes encoding beta-lactamase, sulfonamide (sulI), and tetracycline (tet40) were most common across all sewage (Additional file 4: Dataset 2). Resistance gene profiles indicated distinct seasonal clustering (Adonis test, P < 0.01) (Fig. 1b) and paired-samples t tests further demonstrated significant seasonal differences within most cities (P < 0.05), except the cities of Shenzhen (SZ), Tianjin (TJ), and Xi’an (XA) (Additional file 5: Table S2).

17The relative abundance of antibitotic resistance genes in summer (1.73 copies per bacterial cell) was significantly higher than that in winter (1.15 copies per bacterial cell) (P < 0.01) (Fig. 1c and Additional file 3: Figure S1b). In terms of abundance, winter samples contained highly abundant bacteria (1.21 × 1012 cells/L) (P < 0.01) and resistance genes (1.79 × 1012 copies/L) (P < 0.01), while summer samples were found to harbor lower bacterial abundance (1.70 × 1011 cells/L) and lower resistance gene concentration (3.27 × 1011 copies/L) (Fig. 1d, Additional file 3: Figure S1c and Figure S1d). Significantly different seasonal abundances were observed for 27 ARG subtypes (Additional file 6: Figure s2a).

18No distinct regional distribution pattern of the antibiotic resistome was observed among the sewage samples from different cities (Additional file 7: Figures S3a and S3b and Fig. 2a, b). Based on the demographic data (Additional file 8: Dataset 3), the total volume of domestic sewage discharge ranged from 0.329 to 7.846 million tons/day across major Chinese cities in 2014. The ARG load in the major Chinese cities was calculated with a range from 5.89 × 1012 to 7.85 × 1014 copies/person/day (Additional file 9: Figure S4). The urban ARG burden in Chinese administrative regions was calculated by multiplying the median value (9.47 × 1013 copies/person/day) of ARG load by the urban population, resulting in a range from 5.40 × 1019 to 6.91 × 1021 copies (Additional file 8: Dataset 3). ArcGIS mapping of antibiotic resistance showed significantly higher ARG burden in the east of China, which was 1–2 orders of magnitude higher than those in the west of China. The antibiotic resistance distribution was distinguished by the “Hu Huanyong line”, which delineates a striking difference in the distribution of China’s population (Fig. 3). A strong spatial dependency was observed in the ARG burdens between geographically nearby provinces with Moran’s I index of 0.173 (variance = 0.0056; z score = 2.709; p value = 0.007). Moran’s I index > 0 indicates spatial autocorrelation and larger values of Moran’s I indicate higher spatial autocorrelation. P value < 0.01 indicates an extremely significant spatial autocorrelation.

19From PCR amplicons spanning the V4 and V5 hypervariable regions of the 16S rRNA gene, 6,174,489 high quality sequences (22,149–147,635 per sample) were clustered into 74,138 OTUs (2551–10,691 for each sample, mean = 5,509) (Additional file 1: Table S1). Higher OTU numbers and higher microbial diversity (Chao 1 index, P < 0.01) (Fig. 4a) were observed in summer sewage. Overall, microbial cohorts closely clustered by sampling time (Adonis test, P < 0.01) (Fig. 4b) and sewage microbiomes between seasons were more heterogeneous than those within either season (t test, P < 0.01) (Fig. 4c). A significant distance-decay effect was also observed—similarity in microbial communities between any two cities decreased with increasing geographic distance (r = −0.364, P < 0.01) (Fig. 4d). Similar to antibiotic resistome, no geographical cluster of either bacterial community or shared bacterial OTUs was also observed (Additional file 7: Figures S3c and S3d and Fig. 2c, d).

20Proteobacteria, Bacteroidetes, Firmicutes, and Fusobacteria were the dominant phyla, accounting for 58.7 to 98.5% of sequences within the 116 samples (Additional file 10: Figure S5a). Bacteroidetes had significantly higher abundance in summer samples than that in winter samples (P < 0.01), while both Firmicutes and Actinobacteria were more abundant in winter than those in summer (Additional file 10: Figure S5b). Significantly different seasonal abundances were also observed in several bacterial classes, for example, Deltaproteobacteria, Bacteroidia, and Clostridia (Additional file 6: Figure S2b). At the genus level, the most abundant 30 genera were mainly classified to Proteobacteria, Bacteroidetes, and Firmicutes, and Bacteroides, Prevotella and Acinetobacter were the top 3 abundant genera (Additional file 11: Figure S6a).

21128 resistance genes were shared in more than 80% of samples, accounting for 95.6% of all the ARGs observed. A set of 31 resistance genes were found in all sewage samples, and this core resistome contributed 57.7% (ranging from 29.4 to 84.3%) to the total ARGs detected (Fig. 5). Among the core resistome, the genes for aminoglycoside, tetracycline, beta-lactam, and MLS resistance were dominant (Additional file 12: Figure S7a). Resistance genes sulI, tet40, and one encoding chloramphenicol acetyltransferase were the most abundant resistance subtypes (Fig. 5). The core resistome clustered by season (Adonis test, P < 0.01) (Additional file 12: Figure S7b), but geographical clustering was not observed (Additional file 7: Figure S3 and Fig. 2b).

22Common 16S rRNA OTUs accounted for 64.6% of the total reads in more than 100 samples. A highly shared prevalence of microbiome among all sewage samples was also observed to include 88 classified OTUs, and accounted for 13.6–67.7% (on average 48.8%) of the total bacterial abundance in each sample (Fig. 6a, b). The core OTUs belonged to 33 dominant genera that were affiliated to 7 phyla (Additional file 12: Figure S7c), of which Proteobacteria, Bacteroidetes, and Firmicutes were the most abundant. At the genus level, the most prevalent OTUs were Bacteroidetes, Prevotella, and Trichococcus (Additional file 11: Figure S6b). Similar to the core resistome profiles, these shared OTUs were separated by season without a geographical distribution pattern (Fig. 2d and Additional file 12: Figure S7d).

23Co-occurrence patterns between antibiotic resistance and bacterial assemblages were explored based on strong (ρ > 0.7) and significant (P < 0.01) correlations (Fig. 7). There were more complex and dense correlations between bacterial communities in winter than those in summer, whereas looser relationships with the antibiotic resistome were observed in summer (Additional file 13: Figures S8). Co-occurrence, using network analysis, is summarized in Additional file 14: Dataset 4. Mantel test indicated that the antibiotic resistome was significantly correlated with bacterial phylogeny (P < 0.01). Exploration of connections between ARGs and bacterial genera showed that most genera from the same phylum had similar antibiotic resistance profiles (Additional file 15: Table S3), but this does not prove that bacterial OTUs were actually hosts of resistance genes. In addition, significant correlation between core resistance genes and core bacterial OTUs was also observed (Procrustes test, M 2 = 0.927, P < 0.001, 9999 permutations).

24A total of 205 genera of gut bacteria were detected in 116 samples with abundance ranging from 6.1 to 59.9% (average 28.6%) in each sample. Both the diversity and the relative abundance of gut bacteria in winter samples were higher than those in summer samples (Additional file 11: Figure S6e). A significant seasonal variation of gut bacterial community structure was also observed (Adonis test, P < 0.01) (Additional file 11: Figure S6f). Geographical profiles varied from city to city, with no obvious regional clustering in either total gut bacteria or shared gut bacteria (Fig. 2e, f). Of the detected human gut microbiota, Proteobacteria, Firmicutes, and Bacteroidetes were most abundant phyla. The most abundant gut bacterial genera in sewage were affiliated with Acinetobacter, Arcobacter, and Paludibacter (Additional file 11: Figure S6c). The shared human gut microbiota, including 32 human gut genera detected in all sewage accounted for 49.6% (ranging from 34.6 to 70%) of the shared OTUs (Fig. 6c).

25The urban sewage resistome represents the emission of antibiotic resistance from the gastrointestinal tracts of citizens into wastewater treatment plants. A large-scale sampling of municipal sewage from 17 major cities across China was performed in this study. The seasonal variation and geographical distribution of the urban sewage antibiotic resistome were characterized. Municipal sewage harbored diverse and abundant resistance genes, conferring resistance to almost all antibiotics, highlighting that municipal sewage could be a major conduit for transferring antibiotic resistance genes into the environment. Significant seasonal differences were observed in the urban sewage antibiotic resistome (P < 0.01). Seasonal temperature changes might have a significant influence on the variation of antibiotic resistance and on the composition of microbial community. The seasonality of microbial communities in lakes, soil, and sludge has been well documented. Temperature and temperature-dependent organic matter load can foster the proliferation of microbial taxa that carry resistance genes, or improve the growth of non-resistant microorganisms. Antibiotic administration could also explain the clear seasonality in shifts within the sewage resistome. Seasonal variation in antibiotic consumption driven by associated seasonality in pathologies exerts selective pressure, leading to selection and subsequent dissemination of antibiotic resistance genes in wastewater. Thus, seasonally variable release of antibiotics, bacteria, and resistance genes into municipal sewage can alter bacterial populations and remodel their resistome. A recent study on sewers further supported the speculated reason, demonstrating a clear seasonal pattern in the relative abundances of resistance genes, and that this coincided with the overall rates of antibiotic prescription.

26Higher concentrations of ARGs were detected in winter sewage, being approximately one order of magnitude higher than those of summer. This finding was supported by a recent study, where increases in ARGs in sewers were always encountered in colder seasons of the year, when the more frequent seasonal epidemic diseases contributed to the therapeutic prescription of antibiotics. In addition, dilution of urban sewages by increasing domestic daily water discharges in summer may be another explanation for lower ARG concentration in summer. Similar observations were found in urban streams, where total bacterial numbers were the highest in winter. Although ARGs absolute concentrations in winter sewage were greater, a significantly higher relative abundance of ARGs in summer sewage was observed. The major reason for such disparity might be that the bacterial density in winter sewage (1.21 × 1012 cells/L) was much greater than that in summer sewage (1.70 × 1011 cells/L). Bacterial biomass in sewage has been quantified using flow cytometry within the range from 1010 to 1012 cells/L, and this is consistent with our results when normalizing for 16S rRNA gene copy number per cell.

27The total output of resistance genes was estimated by considering the amount of domestic sewage and urban populations in Chinese administrative districts to quantify the regional ARG burden at a national scale. A strong spatial dependency in the distribution of ARG abundance in various administrative areas were observed, with two main regions separated by the demographic “Hu Huanyong line”, which is based on climatic zonation and population density. It was previously reported that industrialization was correlated with the antibiotic resistance burden of the human gut, and this was in turn driven by age, diet, cultural tradition, climate, pathogen carriage, and periodic perturbation, for example, by antibiotic exposure. A similar geographic distribution was found in the antibiotic emission densities in Chinese river basins, suggesting that human activities are the major driver of resistance gene distribution.

28Geographical clustering was not observed in the structure of either antibiotic resistome or bacterial community in Chinese sewage. The core antibiotic resistome and the core microbial community were stable across WWTPs. The resistome closely correlates with host-related bacterial phylogeny in sewage, indicating that the shared resistome and the core microbiota might play a vital role in the profile of urban resistome and its microbial community. Despite no distinct geographical grouping, there was a distance-decay effect in the similarity in bacterial community composition. This has been specifically reported in freshwater bacterial communities, phyllosphere bacteria, and more generally.

29Median fecal dry mass production is estimated at 29 g per person per day; and human intestinal contents range from 1010 to 1011 bacterial cells per gram (dry weight). Therefore, it was estimated that approximately 1011 ~ 1012 bacterial cells per person per day were discharged into sewage. Given the proportion of these cells that carry antibiotic resistance, there are clear pathways for dissemination of resistance genes via sewage. Although the core resistome was shared by all populations investigated here, there were differences in the abundance of ARGs between urban areas. This suggests that monitoring sewage systems for ARGs could provide a real-time estimate of antibiotic resistance threats in specific areas, and this in turn could be used to inform treatments and to promote stewardship of antibiotics.

30Currently, sound and necessary data on seasonal and geographical characterization of antibiotic resistome in urban sewage is still lacking. This study provided solid evidence for seasonal and geographical patterns of the profiles of antibiotic resistome and potential ARG hosts via a national-scale survey. Seasonal variation in both antibiotic resistomes and bacterial communities was observed in urban sewage. No distinct geographical cluster was found in the distribution of the resistance genes and bacterial community composition. The demographic “Hu Huanyong line” separated the regional ARG burden into two main regions, suggesting human activities might be the major driver of antibiotic resistance burden distribution. A core, shared antibiotic resistome accounted for more than 50% of the total resistance genes, and was significantly associated with the core microbial community. The shared resistome and the shared bacterial community exhibited a distinct seasonal distribution, but did not show geographical clusters, indicating that the share resistome and the core microbiota might play a vital role in the profile of urban resistome and its microbial community. In addition, the strong correlations between resistome and bacterial communities, especially between the core, shared resistome and the core human gut microbiota, indicated the contribution of human gut microbiota to the dissemination of antibiotic resistance. These data provide dynamic background (seasonal and geographical variation) for mitigation activities in WWTPs based on the presence of ARGs and practical guide for improving antibiotic management in the urban sewage.