Bacterial proportion in WWTP effluents
Surface water from the sewage effluent samples was collected from eight WWTP sites along the Tama River and around Tokyo Bay (Fig. 1) and the recreational beach area (BEC1) (Fig. 1). Basic information on the WWTPs is summarized in Table 1. The sampling was conducted in summer and winter seasons for two years between August 2017 and February 2019 (Table S1). The collected water samples were filtered through a 0.2-µm pore membrane, and DNA was extracted from the membrane, followed by metagenomic DNA-seq short-read sequencing (Fig. S1).
Taxonomic classification of domain rank (Fig. S2A) suggested that the effluents of fresh water area (WP1–8) were mostly composed of bacterial sequences, whereas those of brackish water (WP9 and BEC1) were rich in eukaryote organisms, such as mussels (Mytilus edulis) and diatoms (Chaetoceros spp.), found in seawater. Further taxonomic classification of family rank (Fig. 1B) suggested that fresh water area (WP1–8) effluents predominantly contained Comamonadaceae, which are aerobic gram-negative β-proteobacteria as reported in sewage active sludge digesters [8], and Mycobacteriaceae, which are rich organisms in sewage active sludge [9]. The effluents of brackish water (WP9 and BEC1) were rich in Rhodobacteriaceae, which mostly habitats aquatic environments.
Non-metric multidimensional scaling (NMDS) showed that plots of the effluents of fresh water area (WP1–8) were closely scattered, but the plots were markedly divided in a season-dependent manner (Fig. 2A). Linear discriminant analysis effect size (LEfSe) analysis revealed that multiple aquatic habitat bacteria were mostly observed in winter, whereas some specific members of the human gut flora (Streptococcus, Escherichia, and Anaerococcus), Cyanobacteria (Dolichospermum), and purple sulfur bacteria (Rheinheimera) increased in summer (Fig. 2B), demonstrating typical seasonal features.
ARG proportion in WWTP effluents
Abovementioned metagenomic DNA-seq short reads were analyzed by Blastn homology search in the ResFinder ARGs database under normalization using reads per kilobase of exon per million mapped reads (RPKM). Prior to the evaluation of ARGs, we assessed whether all metagenomic DNA-Seq short reads were properly obtained to characterize the bacterial population size using RPKM normalization for 16S rRNA genes and rpoB RNA polymerase β-subunit orthologs (Fig. S2B). All samples, except for BEC1, in winter 2019 showed no significant difference with an average of four copies of the 16S rRNA gene per one copy of rpoB, suggesting that RPKM estimation is an effective approach to characterize the relative abundance of ARGs. Resistance genes to quaternary ammonium compounds (QAC), sulfonamide (SA), aminoglycoside (AMG), β-lactam, and macrolide (MAC) were markedly detected in effluents from all WWTPs, except for WP9 and BEC1 (Fig. 3A). Unlike the bacterial proportion shown in Fig. 2A, a seasonal difference was not observed in ARG variations (Fig. 3B), but it appeared to be dependent on sampling year and season.
We speculated that possible occasional trends of ARB dissemination may affect ARG variations at every sampling and WP site. The detected ARGs in each AMR category were further classified into orthologous genes (or gene families) (Fig. 4A and Table S2). The most detected ARG was the SA resistance gene sul1 at 440.5 RPKM (Table S2); sul2 was detected at rather low levels, but sul3 was not detected. The second most detected ARG was the QAC resistance gene qacE at 203.8 RPKM (Table S2), suggesting that basic gene sets (sul1 and qacE) in the class 1 integron [10] were the most predominant in the detected ARGs.
Regarding clinically important ARGs, multiple AMG- and MAC-resistant genes were identified (Fig. 4A). The seven most predominant AMG resistance genes belonged to the aadA and aph gene families (Table S2), which contribute to streptomycin resistance, implying that streptomycin could be the most potent selective agent among aminoglycosides. The two most predominant MAC resistance genes were msr(E) and mph(E), which are primarily located on the chromosomes of Acinetobacter and Proteus species (The Comprehensive Antibiotic Resistance Database, CARD https://card.mcmaster.ca/ontology/39685).
Regarding β-lactamase gene type, blaOXA, blaGES, and blaIMP were the most markedly detected (Fig. 4A), whereas the blaCTX-M cluster, which is a popular extended-spectrum β-lactamase (ESBL) gene in clinical settings, was in observed in minute abundance at 0.82 RPKM (Table S2), indicating that other β-lactamase type genes, such as blaOXA and blaGES, may be present in baseline levels were derived from environmental bacteria but not pathogenic bacteria.
In addition to the distribution of ARG, metagenomic analysis suggested a unique season-dependent ARG detection. For instance, qnrS2, qnrS6, and aac(6’)-Ib increased markedly in summer (Fig. 4B), whereas aadA13 was detected in winter (Fig. 4C). These season-dependent ARG detection patterns suggested that, to some extent, ARGs could be affected by seasonal elements in a particular temperature-dependent manner (above 30°C in summer and below 10°C in winter; Table S1) similar to that observed in the growth of ARB in activated sludge of the anaerobic–anoxic–oxic (A2O) water treatment system.
However, another possibility remains that seasonal (or temporal) infectious trends may influence the differences in antimicrobial use in the WP-dependent treated area, leading to differential human gut flora as excreta. Thus, possible WP-dependent ARG differences were elucidated using a fold change (log2 ratio) of RPKM (Fig. S3A). In fact, season-dependent ARGs as well as steadily increasing ARGs were observed in WP-specific features (Fig. S3B). WP2 effluents showed that the RPKM of major ARGs [AMG, aadA2 and aadA4; MAC, msr(E) and mph(E); SA, sul1; and QAC, qacF] increased steadily (Fig. S3B), suggesting that the population of major ARB species may increase predominantly. Likewise, the WP-specific increase in ARGs was identified at other WPs, and metagenomic RPKM value could be a suitable index as one of the evaluation standards to monitor such increasing trend of ARB in WWTP-specific areas.
Heavy-metal resistance genes (HMRGs) in WWTP effluents
Resistance genes to mercury (Hg), copper (Cu), and arsenate (As) were markedly detected in effluents from all WWTPs (Fig. S4A and Table S3). Unlike the bacterial proportion shown in Fig. 2A, seasonal difference was not observed in heavy-metal resistance (Fig. S4B), but it appeared to be dependent on sampling year and season as well as on ARGs (Fig. S5).
Genome analysis of ESBL/carbapenemase-producing organisms (EPO and CPO) isolated from WWTP effluents and environment
The metagenomic approach is good for revealing overall ARGs, but it is insufficient to characterize potentially pathogenic ARB. Therefore, we isolated EPO and CPO from fresh 50-mL samples of each effluent (Fig. S6). CHROMagar ESBL selection demonstrated that potential clinically pathogenic Enterobacteriaceae (E. coli, Klebsiella, and Enterobacter spp.) were observed markedly at approximately 500 colonies from the effluents of WP4, WP5, and WP8 at all tested sampling times (Fig. S6). More than hundred typical chromogenic colonies, which would be sufficient to characterize environmental pollution due to EPO and CPO, were isolated at each sampling time (Table 2). These isolates were subjected to whole-genome sequencing (Fig. S6 and Table S4). CTX-M-positive organisms were isolated as follows: 134 (31%) (Table 2), CTX-M-9 group (blaCTX-M-9, blaCTX-M-13, blaCTX-M-14, blaCTX-M-24, blaCTX-M-27, and blaCTX-M-65) was the predominant group with 71 isolates (53%), and the second was the CTX-M-1 group (blaCTX-M-3, blaCTX-M-15, and blaCTX-M-55) with 48 isolates (36%) (Table S4). The EPO and CPO isolates on CHROMagar ESBL and its sequence type (ST) are summarized in Table S5. Eighteen CPO isolates (seven blaIMP-positive, eight blaKPC-positive, and three blaNDM-positive) were identified (Table 2).
Among the 71 E. coli isolates, nine isolates of predominant ST131 (11 isolates) possessed blaCTX-M-27 gene (Table S4). Core-genome single nucleotide variation (SNV) phylogeny (see detailed procedures in Fig. S7) of the ST131 isolates were compared with publicly available E. coli ST131 draft or complete genome sequences (a total of 315 strains; see Table S6), suggesting that primarily four subclonal types of ST131 were identified from the WWTP effluents, two isolates of ST131-fimH41 with blaCTX-M-27, seven isolates of ST131-fimH30 with blaCTX-M-27, one isolate of ST131-fimH30 with blaCTX-M-3, and one isolate of ST131-fimH30 with blaCTX-M-15 (Fig. S8).
One isolate (WP5-S18-ESBL-09) of O16:H5-ST131-fimH41 with blaCTX-M-27 exhibited marked clonality, showing less than 12 SNVs with clinical isolates from several countries (Southeastern Asia, EU, and Oceania) (Fig. 5A), and genome recombination and pan-genome analysis suggested that the genome structure of these strains was very stable, except for some ARGs that were transferred with conjugative plasmid (Fig. 5B). Three isolates of ST131-fimH30 with blaCTX-M-27 exhibited marked clonality, showing 8–13 SNVs to clinical isolate SAMD00126441 in Japan (Fig. S8), suggesting that a clinically potential subclone (O25:H4-ST131-fimH30 with blaCTX-M-27) [11] has been detected in WWTP effluent from Japan. In addition to E. coli ST131, ST38 (10 isolates), ST10 (five isolates), ST405 (four isolates), ST69 (three isolates), and ST648 (three isolates) were observed in this study (Tables S4 and S5). Core genome SNV phylogeny for each ST isolate, including publicly available strains (Table S6), suggested that CTX-M-14 is one of the major CTX-M variants in E. coli ST38, ST69, ST405, and ST648 (Fig. 6). For instance, two isolates of ST69 (WP7-S17-ESBL-01 and WP8-S17-ESBL-12) harbored blaCTX-M-14, but relative clones in other countries harbored blaCTX-M-15, blaCTX-M-24, and blaCTX-M-27 (Fig. 6).
Among the hundred isolates of Aeromonas species (Table 2), 10 isolates of A. hydrophila, 10 isolates of A. caviae, one isolate of A. media, and four isolates of A. veronii were identified as CTX-M producers (Table S4). Although four isolates of CTX-M-14 producing A. caviae were detected in different WWTP effluents (WP2, 4, and 7), core-genome SNV phylogeny revealed that they were clonal in 21 SNVs with same recombination and pangenome regions (Fig. S9A), suggesting that the A. caviae clone may be successfully predominating in the general WWTP environment in Tokyo but not in WWTPs at specific locations.
Among the 51 isolates of Klebsiella species (Table 2), eight isolates of K. pneumoniae, seven isolates of K. quasipneumoniae, two isolates of K. variicola, and seven isolates of other Klebsiella sp. were identified as CTX-M producers (Table S4). Core-genome SNV phylogeny revealed that four isolates of CTX-M-3-producing K. quasipneumoniae ST668 exhibited clonality in 14 SNVs and shared similar recombination and pangenome regions (Fig. S9B), although these isolates were obtained from the same place (WP5) but at different sampling times (summer 2017, summer 2018, and winter 2019), suggesting that the ST668 clone may remain in active sludge at WP5 for at least one and a half years.
Plasmidome analysis of EPO and CPO isolated from WWTP effluents and environment
ESBL and carbapenemase β-lactamase genes are acquired by conjugative plasmid transfer among variable strains in Proteobacteria. Genetic network analysis among the AMR plasmids could be useful to uncover the mode of relational transfer and to trace dissemination. Of the potential 433 EPO/CPO isolates (Table 2), 85 strains were identified as complete genome sequences in this study (Table S4). Of the 304 complete plasmid sequences, 73 β-lactamase-positive plasmids (Table S7) were subjected to plasmidome network analysis based on orthologous analysis (see details in Fig. S10). To characterize global ARG dissemination through plasmid transfer, 758 publicly available complete plasmid sequences showing clear description of isolation and source were selected from 19,904 complete plasmid sequences (from NCBI in November 2019). A total of 831 complete plasmids (73 in this study and 758 in public domain; see Table S8) were subjected to pangenome analysis using Roary and NMDS (vegan in R package) and visualized; the results revealed that 14 communities were classified and clustered (Fig. 7).
Most network communities (Co) showed notable incompatibility (Inc) replicon-dependent distribution, suggesting that the network analysis was well-performed based on the genetic features of each Inc replicon. The major Co comprised Co10, Co11, and Co12 of the shared 423 plasmids among all tested 831 plasmids. Co10 was primarily IncFIA and IncFIB (AP001918) replicon plasmids harboring blaCTX-M gene in E. coli from human sources, whereas Co12 was primarily IncFIA(HI1) and Inc FIB(K) replicon plasmids harboring blaCTX-M, blaKPC, and blaGES genes in Klebsiella from human sources (Fig. 7). This network analysis highlighted the potential mutual role of Co11 plasmids between Co10 and Co12 plasmids because Co11 plasmids are placed around Co10 (rich in E. coli) and next to Co12 (rich in Klebsiella). Co11 is primarily rich in IncFII, Inc FII(pHN7A8), IncX1, and IncN replicon plasmids harboring blaCTX-M, blaKPC, and blaNDM genes and mainly detected in E. coli from human and animal sources. Therefore, Co11 plasmids may contribute to the horizontal ARG transfer between E. coli and Klebsiella. Regarding other notable Co-exhibiting specific replicon and host bacteria, Co3 comprised IncHI2 and IncHI2A replicon plasmids harboring blaCTX-M and blaIMP genes and was mainly detected in Enterobacter from human sources. Co4 comprised IncB/O/K/Z and IncI1-gamma replicon plasmids harboring blaCTX-M gene and were primarily detected in E. coli from human, environmental, and animal sources. From the viewpoint of carbapenemase genes, gene-specific distribution was found as follows: blaIMP in Co1 and Co3; blaKPC in Co6, Co8, and Co12; blaNDM in Co5, Co11, and Co13. Such an ARG distribution could be involved in plasmid replicon and its bacterial hosts.