Comparing against a custom AS reference genome, ARG, and MGE catalogue
We hypothesized that hospital sewage would introduce substantial ARG and MGE microdiversity into the SBRs. To assess this, we developed two genomic catalogues comprising (1) metagenome-assembled genomes (MAGs)/whole genome sequences of isolates recovered from the SBRs or influent and (2) metagenomic scaffolds bearing ARGs or mobileOGs (i.e., MGE hallmark genes) recovered using multiple hybrid assembly, coassembly, and binning strategies (see Methods). The combined metagenomic assembly yielded over 300 Gbp of hybrid assembly data and after binning produced 876 species-level dereplicated medium- or high-quality MAGs (Table S2, Fig. S4).
For subsequent analysis of microdiversity and HGT, we relied on the ARG and MGE context catalogue derived from the 300 Gbp of hybrid assembly data (See Methods) and subsequently related those findings back to potential hosts via analysis of the MAGs. The final resistance gene and mobileOG catalogue consisted of 1,354,363 contigs (total assembly size 1,124,994,800 bp; N50 = 1,587 bp). From these assemblies, a total of 910 unique reference resistance genes were detected (535 ARGs and 375 metal resistance/biocide resistance genes), which is comparable to previous studies.(5)
A diverse array of mobile resistance genes were detected in the MGE and ARG catalogue. Mobile resistance genes were classified according to their co-occurring mobileOGs. The classification scheme was such that individual resistance genes could belong to one or multiple categories of integrative element (IGE), transposable element (TE), plasmid, phage, or conjugative element (CE). A total of 408 unique mobile resistance genes were detected across 1,544 contigs (Fig. S5). Of these, 9 (2%) were co-localized with phage hallmarks, 54 (13%) with TEs; 97 (24%) with IGEs; 108 (26%) with CEs, and 175 (42%) with plasmids. In addition, 79 (19%) were co-localized on contigs with markers for TEs, IGEs, CEs, and plasmids. This latter category of highly mobile RGs included OXA-205, mer and qac genes, mphE, mphA, and aadA, among others (Fig. S5).
Whereas short read profiles suggested a modest difference between hospital- and municipal-sewage, analysis of metagenomic microdiversity via the Kairos assess workflow revealed 12,675 contigs unique to the hospital sewage-blended feed (Figs. 2,3). The Kairos assess branch of the nextflow pipeline leverages a metagenomic assembly catalogue such as the one created here for resistance genes and mobileOGs that spans multiple gene contexts. This workflow was applied to identify distinguishing regions between gene contexts, extract these regions, and query metagenomic reads against them to determine whether a given contig can be found in a sample. Only 260 contigs were exclusively detected in the background municipal sewage, a strong contrast to the hospital effluent (Fig. S6). Of the 12,675 contigs unique to the hospital sewage-blended feed, 232 encoded resistance genes (183 unique), including 31 ARGs and 52 biocide or metal resistance genes that were not detected in 0% hospital sewage influent. For example, the unique contribution of the hospital sewage included 14 contigs encoding macrolide resistance genes msrE, mphE, and tet(39) in a transposon-like setting (Figs. 3C,D); 39 sul1-bearing contigs (Fig. S7); 38 mphA-bearing contigs (Fig. S8) among others (Table S3). The patterns of abundances predicted by short read to contig alignment concurred with the trends observed using Kairos assess (Fig. S8). The overall structure of the mobile resistance gene graphs (Figs. 3A,B) for hospital and native sewage were similar (Fig. S9A,B), except that the hospital-associated graphs spanned fewer distinct taxa relative to the background municipal sewage (Fig. S9C).
AS attenuated the hospital- and background municipal-sewage associated resistome, despite microdiversity level differences in resistance gene fate
The use of the Kairos assess workflow further enabled partitioning of ARG-bearing contigs into hospital sewage-associated and background municipal sewage-associated fractions (Figs. 4, S11). This partitioning allowed us to trace the fate of specific resistance gene contexts by assessing which hospital or native sewage-associated contigs remained detectable in AS after several days of operation. Both the municipal and hospital sewage-associated resistomes were largely attenuated (Fig. 4) when evaluated using either Kairos (Fig. 4A,B) or reads mapping to the contigs directly (Fig. 4C-E). However, different drug classes, as well as individual ARGs within the drug classes, had correspondingly different fates in AS (Fig. 4G). Specifically, there was microdiversity-level variability in the persistence of mphA, msrE, tet(39), tet(G), tet(O), and sul1/sul2 and mer family resistance genes, among others (Fig. S12).
Interestingly, the heterogeneity in gene fate observed here coincided with changes in levels of several antibiotics (Fig. 4H-L). Antibiotics that increased over the course of the experiment included erythromycin (Fig. 4H-K) (Kruskal-Wallis: p = 0.001; Wilcox: median peak area 0 vs. 7.5 ×104, p = 0.005), erythromycin-anhydrous (Kruskal-Wallis: p = 0.009; Wilcox: peak area 2×104 vs. 3.5×104 p = 0.013), and sulfamethoxazole (Kruskal-Wallis: ns; Wilcox: 2×104 vs. 0.9×104, p = 0.04). Increases in erythromycin were concomitant with a near 40-fold reduction in median levels of clarithromycin, another macrolide antibiotic, although this was not statistically significant (Fig. 4H). We note that these antibiotic concentrations were not elevated as a result of the addition of hospital effluent, as the observed changes occurred in all reactors. Instead, the increases reflect changes in the local sewage used as the background feed. Macrolide prescription rates in outpatient settings have been shown to have winter/early-spring peaks(35) suggesting this may have been a byproduct of community antibiotic usage. Therefore, we next addressed the potential consequences of changing antibiotic levels on (1) shifts in the abundance of bacteria harboring the resistance genes, and (2), the potential for HGT of contaminant-associated resistance genes.
Putative hosts of sulfamethoxazole and erythromycin resistance genes display disparate trajectories in the presence of elevated antibiotic levels
If the fate of specific resistance gene contexts were chiefly driven by the persistence of their respective hosts, it would be expected that the relative abundance of the hosts would follow those of the resistance gene contigs. While examination of bin abundances revealed a pattern of change over the duration of the experiment that mirrored fluctuating antibiotic contaminant levels (Fig. 5A,B), there was no evidence of bulk selection by erythromycin or sulfamethoxazole for hosts of these resistance genes specifically (Fig. 5C). However, different putative hosts displayed disparate trajectories (Fig. 5C-E). For example, bin39 (phylum Myxococcota and genus Nannocystis) displayed changes in relative abundance consistent with an enriching effect (Wilcox: median 500 RPKM vs. 1,000 RPKM, p < 0.001) (Fig. 5D) and bin96 (phylum Proteobacteria and genus Thauera) displayed a moderate reduction (Wilcox: median 75 RPKM vs. median 45 RPKM, p < 0.001) followed by an increase in relative abundance (Wilcox: median 45 RPKM vs. 75 RPKM, p < 0.001) (Fig. 5E). By contrast, bins 115, 107, and 86 decreased or remained unchanged over the duration of the experiment (Fig. 5G). Whereas multiple bins bearing macrolide resistance genes were found, only one bin, bin82 (phylum Proteobacteria and genus Dokdonella) bearing a sulfonamide resistance gene (sul2) was found. The sole MAG with sul2 displayed decreased relative abundance over time (Wilcox: median 450 RPKM vs. 200 RPKM, p < 0.001). Further, relative abundance profiles of persisting or attenuated contigs (as predicted by Kairos) followed expectations in that attenuated contigs corresponded to hosts that decreased in AS relative to influent and persisting contigs corresponded to hosts with abundances that remained the same or increased (Figs. S13, S14).
Multiple potential pathways underly putative in situ transfer of sulfonamide and macrolide resistance genes in the presence of elevated antibiotics
It was observed that resistance genes were found across diverse genera, classes, and phyla (Fig. 2A-C), suggesting the potential for gene sharing and possibly HGT. Because of the changes in profiles of antibiotics, we assessed the potential for in situ HGT possibly linked to the antibiotics using the general framework proposed previously,(25) with formal and case-specific hypotheses crafted for this experimental design (Fig. 6A,B). In this case, in situ HGT strictly refers to any occurrence of cross-taxa gene sharing with a pattern of presence/absence in samples consistent with an HGT event, or an enrichment of a preexisting HGT, in the sampled period of time.
For this analysis, we applied the Kairos derep-detect workflow to identify contigs for which identical resistance gene or mobileOGs were found, but where different taxonomic assignments were predicted (i.e., gene sharing or potential instances of HGTs). This suggested the potential for extensive HGT across multiple taxonomic levels including phylum (n = 4,919), (Figs. 6C, S15, Table S4), class (n = 1,884) (Table S5), order (n = 3,488) (Table S6), family (n = 983) (Table S7), and genus (n = 1,400) (Table S8). Genes shared between phyla included resistance genes (n = 143) APH(6)-Id, APH(3'')-Ib, OXA-205, qacH, ermG, ermB, mel, mphE, msrE, mphA, mphF, sul2, tet(C), and tet(39) (Fig. 4B). The majority (n = 3,413) corresponded to mobileOGs of diverse categories (Fig. 6C). Gene sharing network analysis revealed dense linkages connecting Proteobacteria to Bacteroidota and Actinobacteriota, but not Actinobacteriota to Bacteroidota (Fig. 6D), These connections were negatively correlated with GC content dissimilarity (Spearman’s rank rho = -0.24, p < 0.001). Of the 12,685 potential HGTs across all taxonomic levels, 1,608 met in situ criteria, including transfers of mphA, sul1, and sul2 (Fig. 6D) and 14 other resistance genes (OXA-2, qacH, sul1, OXA-205, merF, merP, merT, qacF, merA, merD, merE, merR2, and merT). We focused subsequent analysis on putative transfers of mphA and sul1/sul2 as these genes were disproportionately persistent (Figs. 4G, S12) and were most likely to inform the potential impact of erythromycin and sulfamethoxazole.
Analysis of predicted recipients and donors suggested a discrete set of potential transfer pathways (Fig. 6E,F). For mphA, a single recipient in the class Polyangia (phylum Myxococcota) was predicted for two potential donors of the orders Xanthomonadales (phylum Proteobacteria) and Myxococcia (Myxococcota) (Figs. 6E, S16). By contrast, the pathway for sul2 included multiple potential recipients and donors and could variably be described by HGT through some route spanning Gammaproteobacteria, Alphaproteobacteria and Bacteroidia. The use of the in situ criteria did rule out transfer from Alphaproteobacteria to Bacteroidia, but not from Gammaproteobacteria to Bacteroidia (Figs. 6F, S17). A single pathway of sul1 transfer was detected that suggested a transfer from Sphingomonadaceae of Alphaproteobacteria to Gammaproteobacteria, however the putative donor and recipients aligned against one another at their respective edges (Fig. S18).
HGT of mphA encoded by a novel myxophage highlights complex ecological interactions shaping resistance gene fate
Further scrutiny of the Myxococcota genetic contexts of mphA revealed that they likely were derived from a phage (Fig. 7). This finding is notable, as the role of phages in the evolution of antibiotic resistance remains unclear, especially in environmental matrices. Further inspection revealed that closely-related phage genomes were detected in many samples, but the assemblies were fragmented into two approximately 20,000 base pair segments (Fig. 7A), essentially dividing the genome into halves, with only a few of the fragments encoding mphA (Fig. 7B). Draft genomes were constructed by scaffolding contigs based on their alignment to a similar prophage region in NCBI (Supplementary Methods S1, Fig. S19) (Fig. 7D). The assemblies produced by HybridSPAdes and OPERA-MS predicting the encoding of mphA in the phage genome were validated by the identification of nanopore reads that aligned to both Myxococcota genomes and Enterobacterales plasmids in NCBI (Extended Data 1, Fig. S19). Similar putative prophage regions were detected in numerous publicly available genomes (Fig. 7C, Table S9), but lacked any known resistance genes. However, one putative prophage detected in Danish AS-derived MAG CP064980.1 (Fig. 7D) was integrated near genes encoding macrolide ATP-binding transporter/permease, MacA/MacB (Fig. S20). Relative abundances of putative donors and recipients, and the taxonomic assignment of contigs, suggested that the Proteobacteria context was possibly associated with bin82, a Dokdonella sp. The contig had remarkably similar abundance profiles and concordant taxonomic assignments (Fig. 7F, G). The potential donor of the predicted class-level HGT of mphA, taxonomic assignment Myxococcia, was represented by multiple contigs (Fig. S16). Only two MAGs of class Myxococcia were found, which displayed opposite trends in abundance over time (Fig. 7H,I). Potential recipients of mphA encoded by the myxophage (according to class-level taxonomic assignment ascribed to the phage contig) included 43 Polyangia dereplicated bins (Table S10) spanning all four CAGs.
While the precise pathway of HGT is uncertain, the fact that one of the Myxococcia MAGs displayed trends in abundance similar to those of the Dokdonella sp. suggests that there may have been an ecological linkage between the two populations represented by the bins (Fig. 7E,F). Additionally, one of the putative recipients displayed a genetic context suggestive of prophage integration (Fig. S21). Members of the phylum Myxococcota span a wide range of different environments, including sewage, soil, and marine environments.(36, 37) Most, but not all, have been demonstrated to have some degree of predation.(38) We examined genomic evidence of a predatory lifestyle (e.g., secretion systems and antibiotic biosynthesis pathways(39)) in the genome of the putative donor, bin204, of the class Myxococcia and species Archangium (Fig. 7G). Functional annotations suggested the presence of partial or complete type 1 secretion systems (T1SS), T2SS, T3SS, T4SS, and T7SS; at least one antibiotic biosynthesis monooxygenase; one polyketide synthase; and 19 separate CAZy classified carbohydrate active enzymes (Table S11).