Antibiotic Contamination and AMR Dynamics in the Urban Sewage Microbiome: Insights from a Longitudinal surveillance

Deepjyoti Paul 1

Daizee Talukdar 1

Ramani Shyam Kapuganti 1

Vaishali Gupta 1

N Lekshmi 1

Pradipta Jana 1

Pawan Kumar 1

Jyoti Singh 1

Shalini Kumari 1

Chandana Basak 1

Kajol Kamboj 1

Susmita Bakshi 1

Shruti Lal 1

Subhash Tanwar 1

Roshan Kumar 1

Prabhakar Babele 1

Manish Bajpai 1

Yaswant Kumar 1

Ankur Mutreja 2

Sukhendu Mandal 3

Nitya Wadhwa 1

Sanjay K Banerjee 4

Bhabatosh Das 5✉ Phone+91-129-2876471 Emailbhabatosh@thsti.res.in

1 Centre for Microbial Research BRIC-Translational Health Science and Technology Institute 121001 Faridabad India

2 Cambridge Institute of Therapeutic Immunology & Infectious Disease (CITIID), Department of Medicine University of Cambridge Cambridge United Kingdom

3 Department of Microbiology University of Calcutta 87/1, College Street 700 073 Kolkata India

4 Department of Biotechnology National Institute of Pharmaceutical Education and Research (NIPER-Guwahati) 781101 Changsari, Guwahati Assam India

5 Functional Genomics Laboratory, Centre for Microbial Research BRIC- Translational Health Science and Technology Institute, NCR Biotech Science Cluster, Faridabad – Gurgaon Expressway 3rd Milestone PO box #04 121001 Faridabad Haryana India

Deepjyoti Paul¹, Daizee Talukdar¹, Ramani Shyam Kapuganti¹, Vaishali Gupta¹, Lekshmi N¹, Pradipta Jana¹, Pawan Kumar¹, Jyoti Singh¹, Shalini Kumari¹, Chandana Basak¹, Kajol Kamboj¹, Susmita Bakshi¹, Shruti Lal¹, Subhash Tanwar¹, Roshan Kumar¹, Prabhakar Babele¹, Manish Bajpai¹, Yaswant Kumar¹, Ankur Mutreja², Sukhendu Mandal³, Nitya Wadhwa¹, Sanjay K Banerjee⁴, Bhabatosh Das^1@

1-Centre for Microbial Research, BRIC-Translational Health Science and Technology Institute, Faridabad 121001, India

2-Cambridge Institute of Therapeutic Immunology & Infectious Disease (CITIID), Department of Medicine, University of Cambridge, Cambridge, United Kingdom

3-Department of Microbiology, University of Calcutta, 87/1, College Street, Kolkata-700 073, India

4-Department of Biotechnology, National Institute of Pharmaceutical Education and Research (NIPER-Guwahati), Changsari, Guwahati, Assam 781101, India

@ Correspondence:

Bhabatosh Das: Functional Genomics Laboratory, Centre for Microbial Research, BRIC-Translational Health Science and Technology Institute, NCR Biotech Science Cluster, 3rd Milestone, Faridabad – Gurgaon Expressway, PO box #04, Faridabad – 121001, Haryana, India. Phone: +91-129-2876471; Fax: 0129-2876500. Email: bhabatosh@thsti.res.in

Abstract

The emergence and spread of antibiotic resistance (AMR) in clinically important bacterial pathogens severely compromised the effectiveness of commonly used antibiotics in healthcare. Acquisition and transmission of AMR genes (ARGs) are often facilitated by sublethal concentrations of antibiotics in a microbially dense environments. In this study, we used sewage samples (n = 371) from six Indian cities between June and December 2023 to assess the concentration of eleven antibiotics, microbial diversity, and ARG richness. Our findings revealed the presence of antibiotics from seven drug classes and over 2000 bacterial ASVs. Metagenomic (n = 220) and isolated genome sequences (n = 305) of aerobic and anaerobic bacterial species identified 82 ARGs associated with 80 mobile genetic elements (MGEs). These MGEs were predominantly found in multidrug-resistant bacterial pathogens. Comparative core genome analysis of multi drug-resistant (MDR) bacterial isolates (n = 7166) showed a strong genetic similarity between sewage-derived strains and clinical pathogens. Our results highlight sewage as a significant reservoir for ARGs, where genetic exchanges occur, facilitating the evolution and spread of AMR pathogens within both community and healthcare settings. Additionally, the dipstick-based assay developed for detection of ARGs in the present study could be employed for sewage surveillance in low resource settings for better understanding of resistance prevalence.

Key Words

Antibiotics

Antibiotic resistance

Mobile genetic elements (MGEs)

Pathogen

Microbiome

Sewage

Antimicrobial resistance (AMR) is a critical and urgent global health crisis, threatening human, animal, and environmental health due to the increasing prevalence of multidrug-resistant (MDR) bacterial pathogens¹. The emergence of MDR microbes is driven by the widespread use of antibiotics across various sectors, which exerts selective pressure that favors resistant variants over susceptible microorganisms^2–3. Rapid advances in DNA sequencing technologies have enabled scientists worldwide to decode the metagenomes and whole genomes of diverse microorganisms, facilitating the monitoring of AMR genes (ARGs) and their genetic associations with mobile genetic elements (MGEs). These elements play a crucial role in the rapid spread of resistance through horizontal gene transfer (HGT) among both closely and distantly related microbial species.

The growing need for comprehensive global genomic pathogen surveillance has become increasingly urgent as pathogens continuously evolve. The recent COVID-19 pandemic, one of the most devastating respiratory infectious diseases, highlighted the critical importance of real-time global surveillance systems to detect and track pathogens. Recently, the World Health Organization (WHO) in close collaboration with the Food and Agriculture Organization (FAO) and World Organization of Animal Health (OIE) has developed the One Health AMR Global Action Plan (GAP), to mitigate the AMR crisis. This tripartite organization (FAO-OIE-WHO) is tasked to oversee the development and execution of National Action Plans (NAPs) and strategies which aim at reversing AMR⁴. The WHO is also dedicated to overseeing water, sanitation, and hygiene (WASH) efforts as well as wastewater management, which are crucial for mitigating and reducing the impact of AMR in sewage. An in-depth study on the AMR patterns and burden in sewage and leveraging this data can significantly advance WASH strategies, enabling precise tracking and early detection of AMR pathogens.

Sewage water, a complex and often poorly regulated environment, harbors diverse microbial communities that facilitate gene transfer, which can ultimately reach humans through various pathways. Recent advancements in molecular techniques have significantly enhanced our understanding of AMR in wastewater, allowing for unprecedented insights into microbial communities and their resistomes. However, critical knowledge gaps remain regarding the factors driving HGT in wastewater and the role of MGEs in the active dissemination of ARGs among bacterial species. In this study, we leverage cutting-edge, multi-faceted technologies to explore the dynamics and mechanisms driving AMR in sewage water, a crucial but under-explored area of research.

Results

Quantification of antibiotic contamination in urban sewage using QTRAP mass spectrometer

During June to December 2023, a total of 371 sewage samples were collected from six different cities located across the country, including the national capital region (NCR) of India. A total of 100 sewage samples, collected across multiple rounds from various cities, were analyzed to quantify the concentrations of eleven antibiotics representing seven distinct drug classes using QTRAP mass spectrometer (Fig. 1). The calibration curves generated using pure drug compounds were linear across the tested concentration ranges from 0.1ng/ml to 500ng/ml. The concentrations of tetracycline, amoxicillin, and spectinomycin in sewage were found to be substantially higher (> 15ng/ml) compared to other antibiotics. We observed the highest fluctuation of concentration in spectinomycin in the sewage among all the eleven antibiotics, depending on the site and date of sample collection. Minimal concentration was observed for ampicillin and polymyxin-B (< 0.5 ng/ml). Precision and accuracy data generated in this study were included in Supplementary Table 1a-k.

Fig. 1

Occurrence of various classes of antibiotics in sewage water

This figure illustrates concentrations of different antibiotic residues in sewage water across multiple sampling sites. Each dot represents the concentration of specific antibiotics belonging to macrolide, rifampin, tetracycline, beta-lactams, polymyxin, aminoglycosides, and quinolone group of drugs. The image was prepared using the paid version of Graph pad Prism version 10. The significance of the antibiotics detected were determined using t-test and marked as non-significant (ns), or highly significant (marked with *).

High-throughput longitudinal monitoring of sewage resistome using dipstick-based molecular assay

Given the widespread presence of multiple antibiotics in sewage across the cities, we expanded our analysis to screen for ARGs (n = 16) using a culture-independent dipstick-based molecular assay developed in our laboratory. Our findings revealed significant heterogeneity in the prevalence of ARGs across different cities (Fig. 2). Notably, three ARGs-bla_OXA, qnrS and mphA-were present in more than 30% of the samples tested in this study. The lowest prevalence was observed for the bla_KPC allele (2%), followed by bla_NDM (3%). Other prominent ARGs detected including bla_TEM (28%), tet(26%), bla_CTX−M (21%), cat (19%), and sul1 (16%) (Supplementary Table 2).

Fig. 2

Detection of antibiotic resistance genes (ARGs) using dipstick assay. Heatmap represents the presence (blue color) or absence (grey color) of the sixteen targeted ARGs. X-axis represents the name of the ARGs and Y-axis represents the selected samples for this assay.

Metagenomic investigation of urban sewage microbiomes diversity

A total of 220 sewage samples were analyzed through metagenomic profiling, and sequencing of the 16S V3-V4 region generated 16,679,255 high-quality reads. The number of reads per sample ranged from a minimum of 25,762 to a maximum of 123,274 after quality filtering. The average read count was 75,815, resulting in 1,29,389 Amplicon Sequence Variants (ASVs). (Supplementary Table 3). In our analysis, we have taken samples from different cities within India viz. Faridabad (n = 172), Guwahati (n = 16), Saharanpur (n = 11), Jharkhand (n = 10), Kolkata (n = 7), and Rishikesh (n = 4). Alpha diversity indices (Observed, Chao1, Shannon, and Simpson) were calculated and compared across the six cities. Faridabad district showed a significantly higher alpha diversity index (p-value < 0.05) compared to the other five cities (Fig. 3a). The bacterial diversity in sewage microbiome in Faridabad district is relatively rich, with observed indices ranging from 156 to 879, Shannon indices ranging from 4.15 to 6.4, and Simpson indices ranging from 0.0419 to 0.0033. Beta diversity analysis revealed significant differences between cities, with Bray-Curtis distance matrices demonstrating variances of 19% (Axis 1) and 15.2% (Axis 2) on the ordination plot (Fig. 3b). Permutational Multivariate Analysis of Variance confirmed significant variance among the cities (R²=0.412, p < 0.01).

Fig. 3

Diversity measures and distribution of sewage microbiome based on amplicon sequencing

a) Box plots showing different alpha-diversity measures between the sewage water sites from six different district of India i) Observed ii) Chao1 iii) Shannon iv) Simpson index; statistical significance was observed between Faridabad and all the five sewage except Rishikesh (Wicoxon, p value < 0.0.5). The box indicates the interquartile range (IQR). The median value is represented as a line within the box and whiskers extend to the extreme value i.e. within 1.5 × IQR. b) Beta diversity represented by principal coordinate analysis plots to visualize the difference in sewage microbiome structure between the six sites of India. Each dot represents a particular sample, which is color‐coded based on the sites of the current study. Statistical significance of variance was calculated using permutational multivariate analysis of variance, which reflected that the sites had a variance among them (R²=0.412, p < 0.01). c) Heatmap showing the distribution of top bacterial taxa in the species level in all the six sites. The color code is according to the relative abundance of the bacterial taxa. d-e) Comparisons of the sewage microbiome between community settings and hospital outlets; d: genus level; e: species level comparison.

The species with the most assigned sequencing fragments include Aliarcobacter cryaerophilus across the six cities, except in Jharkhand, where Aeromonas caviae was more prevalent. Several other bacteria with known pathogenic potential, including Escherichia coli, Pseudomonas aeruginosa, Enterococcus species, and Acinetobacter baumannii, were also found to be prevalent in sewage samples collected from various cities across India (Fig. 3c). We compared the sewage microbiome profiles between community settings and hospital outlets in Faridabad district using longitudinal samples collected over seven months. The study included five community sites Badkhal village (n = 27), Ballabgarh (n = 27), Dabua colony (n = 27), Parvatiya colony (n = 26), and Sanjay colony (n = 26), and two hospital outlets: BK hospital (n = 21) and ESIC hospital (n = 18). At the genus level, Aliarcobacter was notably elevated at both hospital outlets, except for ESIC hospital (6.3%), where its relative abundance was lower compared to BK hospital (20.5%). In hospital outlets, the relative abundance of genera such as Aeromonas, Escherichia and Trichococcus increased, while Rombustia decreased compared to community settings. Other Gram-negative pathogens like Klebsiella, Acinetobacter, and Pseudomonas were also prevalent in both environments (Fig. 3d). Specifically, at the species level, Trichococcus flocculiformis and Aeromonas veronii were more abundant in hospital sewage, while the abundance of A. cryaerophilus remained similar across both settings, with a slight decline at the ESIC outlet.(Fig. 3e).

In our shotgun analysis, we sequenced a subset of samples from community settings (n = 37) and hospital outlets (n = 11) to study the taxonomy and diversity indices. Significant differences (p-value < 0.001) in alpha diversity indices (Shannon, Chao1, and Simpson) were observed between the two settings (Supplementary Fig. 1a, i-iii). Aliarcobacter cryaerophilus was the most dominant taxa in sewage water across both sequencing strategies. Taxonomic analysis revealed an enrichment of A. veronii in sewage outlets of hospitals, consistent with amplicon sequencing data (Fig. 4a). Additionally, pathogenic bacteria such as E. coli, K. pneumoniae, K.variicola, P. aeruginosa, A. baumanii, and Enterobacter asburiae were prevalent in the sewage wastewater of Faridabad district in both sequencing approaches (Fig. 4b).

Fig. 4

Shotgun based metagenomic and sewage resistome profiles from community and hospital sites

The heatmap illustrates the distribution of sewage microbial taxa at a) the genus level and b) the species level, with the color code representing the relative abundance of the bacterial taxa. Additionally, there is a heatmap displaying the distribution of sewage resistome data in the community and hospital sewage outlets for c) antibiotic resistance genes (ARGs) and d) efflux genes. e) Network analysis of bacterial species and their associated ARGs is also provided. The red rectangular boxes indicate bacterial taxa, while the differently colored ovals represent the antibiotic drug classes associated with the ARGs.

Seasonal and temporal variations in sewage microbiomes

Monthly sampling from June to December 2023 revealed seasonal shifts in the sewage microbiome. The clustering of points within distinct eclipses in the PCoA plot of each month indicates temporal variation in microbial community composition (R2 = 0.2127, p-value = 0.001). Samples from September had very different scores on the axes and are dissimilar in their microbiome composition with PCoA1 scores greater than 0. The microbiome composition in the months of June and July cluster together, indicating similar community structure during this period. Typically, samples from November and December harbored the most distinct community compositions as evidenced by their separation along the PCoA axes (Fig. 5a). Alpha diversity (Observed richness, Shannon, and Simpson index) was calculated to get a comprehensive view of microbial diversity in sewage samples across months and seasons. In the case of Observed index, the highest richness is observed during monsoon months (June, July, August), with a decline starting from September (post-monsoon) and continuing into the winter months (November & December) (Fig. 5b).The Shannon index shows a similar trend to observed richness but captures both richness and evenness. It peaks during monsoon (June to August) and declines sharply in post-monsoon (September) and further into autumn and winter (November and December) (Fig. 5c). Monsoon months show significantly higher evenness, as no single species dominates (p-value < 0.05).The Simpson index remains relatively stable during the monsoon season but shows significant decreases starting in the post-monsoon season and through winter. This indicates an increase in species dominance as seasons progress, with fewer species maintaining significant populations (p-value < 0.05) (Fig. 5d).

Fig. 5

Changes in the diversity and composition of the sewage microbiome across different months (From June to December 2023)

a) A PCoA (Principal Coordinate Analysis) plot was generated to visualize variations in the sewage microbiome structure across different months. Each dot corresponds to a specific sample, with colors representing the respective months. The statistical significance of variance was assessed using permutational multivariate analysis of variance (PERMANOVA), revealing a significant variance among the sites (R2 = 0.2127, p-value = 0.001). Box plots illustrating different alpha-diversity metrics include b) Observed, c) Shannon, and d) Simpson indices, with statistical significance observed. The box represents the interquartile range (IQR), the median value is indicated by a line within the box, and whiskers extend to the extreme values (within 1.5 × IQR). e) A heatmap depicts the distribution of highly prevalent bacterial taxa at the species level across different months (from June to December 2023). The color code corresponds to the relative abundance of the bacterial taxa

We also examined how the abundance of bacterial taxa varied with the changing seasons. The relative abundance of A. veronii and T. flocculiformis significantly decreased during the monsoon, followed by a gradual increase in abundance as the pre-winter season approached, reaching higher levels during the winter. In contrast, Bifidobacterium pseudocatenulatum, Streptococcus equinus, and Bifidobacterium adolescentis were less abundant during the peak monsoon but increased in abundance during post-monsoon, before experiencing a subsequent decline. Aliarcobacter ellisi, Aliarcobacter faecis, P. aeruginosa, Tolumonas auensis, and Thauera selenatis were abundant during the monsoon and continued to rise through post-monsoon, though their levels dropped during the pre-winter season. Notably, Pseudomonas lundensis was completely absent during both the monsoon and post-monsoon seasons, however it showed increased abundance during the winter season (Fig. 5e).

The resistomes diversity and dynamics in urban sewage

We aimed to analyze the resistome profile in both the community and hospital outlets sewage water sites by using assembled metagenomic sequences. On analyzing, we obtained a total of 629 bins derived from 48 samples used for shotgun sequencing and of these, 244 bins were identified to harbor ARGs. We identified a total of 94 distinct ARGs encoding resistance to 16 drug classes having enzyme mediated or target alteration mechanisms and 66 efflux genes, with varying abundances across the two different sites. The majority of the genes belonged to the beta-lactamase, aminoglycoside, and tetracycline classes. Notably, among all the 94 ARGs, mphE, msrE, mphA, and mphB from the macrolide class were the most common across the samples, while bla_OXA−427, ampC, ampC1, ampH, and bla_MOX−6 were prominent from the beta-lactamase class (Supplementary Table 4a). Tetracycline class gene tetM was also widely detected. Additionally, qnrS1 and catQ from the fluoroquinolone and chloramphenicol classes were observed in multiple samples across community and hospital sites. mphE, msrE, qnrS1, and mphA were found at higher levels in community sewage sites, while bla_OXA−427, dfrG, IsaE, ermT, qnrB4, bla_OXA−12, and tetQ were more prevalent in hospital settings. (Fig. 4c). We identified 66 efflux genes, with CRP being the most common across all samples, while the mexF gene was more prevalent in hospital settings (Supplementary Table 4b, Fig. 4d).

Network analyses of bacterial species and their associated ARGs frequently uncover co-occurrence patterns, suggesting which bacterial species are likely to contain or share resistance genes. A network graph depicting bacterial species-ARG association showed clear separation according to species-level taxonomy (Fig. 4e). Proteobacterial species such as K. pneumoniae, E. coli, P. aeruginosa, A. baumannii, Enterococcus faecium, and Enterobacter spp. harbored a variety of ARGs. Some ARGs like aac(3)-IId, aph(3’’)-Ib, aph(6)-Id, mcr-9, qnrS1, bla_TEM, mphA, mphE, msrE and sul2 were shared among them, while others were uniquely seen in their respective species-specific. E. coli, in particular, had the highest number of assigned ARGs and shared many with other Gram-negative bacteria including K. pneumoniae, Proteus mirabilis, A. baumannii, P. aeruginosa, Enterobacter spp and others (Fig. 4e). We also identified mcr-9, a peptide antibiotic, within the contigs of several Gram-negative pathogens, including E. coli, Enterobacter mori and E. cloacae. Among Gram-positive bacteria, Enterococcus had the most diverse ARGs, which were shared with other Gram-positive species like Staphylococcus aureus, Bacillus paranthracis, and several Streptococcus species (Fig. 4e)

Functional metagenomic analysis of the sewage samples has revealed several pathways likely involved in the degradation of xenobiotic pollutants. Using Illumina shotgun sequencing, we examined the biodegradation-related KEGG pathways in community and hospital sewage sites in Faridabad. Wherein, we identified 17 KEGG categories associated with various xenobiotic pollutants, including chlorocyclohexane, chlorobenzene, benzoate, dioxin, xylene, toluene, polycyclic aromatic hydrocarbons, chloroalkane/chloroalkene, naphthalene, aminobenzoate, nitrotoluene, and steroid, with benzoate being the most prevalent and significant pollutant detected (Supplementary Fig. 1b).

Antimicrobial resistance phenotypes of aerobic and anaerobic bacterial taxa isolated from urban sewage

Culturomics analysis of the samples identified a total of 964 bacterial isolates, encompassing both Gram-negative and Gram-positive from sewage water. The bacterial isolates include MDR strains with aerobic and facultative anaerobic traits. Among the Gram-negative bacteria, notable isolates include E. coli (n = 264), Morganella morganii (n = 127), Proteus mirabilis (n = 72), K. pneumoniae (n = 70), Enterobacter spp. (n = 22), Pseudomonas spp. (n = 18), and A. baumannii (n = 8). In contrast, the Gram-positive bacteria identified are mostly facultative anaerobes, with E. faecium being the most prevalent (n = 138) (Supplementary Table 5a).

All the aerobic isolates (n = 627) and facultative anaerobic isolates (n = 54) included in WGS analysis were subjected to susceptibility testing against 32 antibiotics from 10 different classes (Supplementary Table 5b, Supplementary Table 6a). A majority of these isolates (93.8%; n = 639/681) demonstrated resistance to more than 10 antibiotics (Fig. 6a). Notably, Gram-negative ESKAPE pathogens especially K. pneumoniae, A. baumannii, P. aeruginosa and Enterobacter spp. exhibited resistance to multiple antimicrobial classes such as aminoglycosides, beta-lactams, fluoroquinolones, macrolides, and tetracycline. Among Gram-positive bacteria, most E. faecium isolates were resistant to most of the tested antibiotics, including erythromycin and vancomycin (Fig. 6b, Supplementary Table 6b).

Fig. 6

Resistant profile of sewage isolates against multiple classes of antibiotics

a) Whisker plot illustrates the sewage isolates resistant to the total number of antibiotics, b) The heatmap represents the overall percentage of resistance among bacterial isolates to various classes of antibiotics tested. The white-blue gradient in the heatmap indicates resistance levels ranging from 0% (susceptible) -100% (resistant).

In accordance with the WHO's 2024 guidelines, the study isolates were categorized into critical and high-priority groups. All the A. baumannii study isolates were classified as critical priority pathogens (CPP), due to their resistance to at least one carbapenem antibiotic. Furthermore, among the tested enterobacterial isolates, 69.3% were found to be carbapenem-resistant which includes Enterobacter spp. (78.8%), Citrobacter spp. (72.7%), Proteus spp. (72.2%), K. pneumoniae (70.5%), Morganella spp. (87.8%) and E. coli (57.7). However, 84.7% of these Enterobacterial isolates exhibited resistance to cephalosporins that are also categorized as CPP. Conversely, the high-priority pathogens identified in this study included vancomycin-resistant E. faecium, which accounted for 18.8% of the tested isolates, and carbapenem-resistant P. aeruginosa, with a concerning resistance rate of 88.8% (Supplementary Table 5c).

Genomic insights into MDR bacterial pathogens isolated from urban sewage

A total of 305 isolates cultured from wastewater and hospital sewage oulets were sequenced and analyzed to understand their genetic relatedness and genetic signatures. These isolates belonged to 16 different genera and 31 different species; predominantly the Gram-negative Escherichia (n = 69), Klebsiella (n = 40), Morganella (n = 33), Providencia (n = 31), Citrobacter (n = 29), Proteus (n = 25), and the Gram-positive Enterococcus (n = 22). The overview of all the isolates is shown in the phylogenetic tree based on the 5S, 16S and 23S rRNA sequences (Fig. 7a & Supplementary Table 7a). Global phylogeographical analysis was done to understand the relatedness of the Indian isolates from the sewage and the global representative collection of clinical isolates including E. coli, K. pneumoniae, A. baumannii, P. aeruginosa, E. faecium, and the emerging pathogens M. morganii, P. mirablis, Citrobacter werkmani, and Providencia rettgeri from 2019–2023 (Supplementary Table 7b-j). The phylogenetic trees were constructed based on the core-gene SNPs (Fig. 7b-f & Supplementary Figs. 2a-d). The phylogeographical analysis of all the 2769 E. coli isolates (2700 clinical and 69 sewage) showed that these isolates were quite diverse (Fig. 7b) and widely distributed into eight different phylogroups (A, B1, B2, C, D, E, F, & G). However, the study isolates from sewage belonged to the phylogroups A, B1, B2, C, D, & F and the most prevalent were the phylogroups B1 (n = 22) and A (n = 20). Although B2 isolates were relatively more common in the global clinical cases, this phylogroup was not as prevalent in the Indian clinical cases or the sewage (n = 3). The sewage isolates belonged to 22 different STs/lineages, the predominant were ST167, ST410, ST448, ST648 and ST5869. A few of these isolates that belonged to ST167, ST2083, ST410 and ST448 clustered with the Indian clinical isolates.

Fig. 7

Phylogenetic trees showing genetic relatedness of the isolates

a) Maximum likelihood phylogenetic tree based on the similarity of the 16S, 23S and 5S rRNA sequences gives an overview of the isolates from sewage. Phylogenomic trees based on the core-gene SNPs shows the genetic relatedness of sewage isolates with the global clinical cases show b) E. coli (n = 2769), c) K. pneumoniae (n = 1314), d) A. baumannii (n = 697), e) P. aeruginosa (n = 780), f) E. faecium (n = 660). The trees are annotated with metadata such as sample source, geographical location and year of collection. The scale bar is indicated in the figures.

The phylogeographical analysis of all the 1314 K. pneumoniae isolates (1281 clinical and 33 sewage) showed these isolates to be less diverse and distributed into distinct clades (Fig. 7c). Although the study isolates are distributed on different clades, these isolates clustered closely with the Indian clinical isolates. The most prevalent STs in sewage isolates were ST15, ST37, and ST16 which were also observed to be widely isolated from clinical cases in the US, East and South-East Asia. Although not prevalent, ST11, ST101, ST231 and ST147 found in sewage isolates were also observed in Indian clinical cases. The phylogeographical analysis of 697 A. baumannii isolates (689 clinical and 08 sewage) showed that the isolates are clonal (Fig. 7d). The sewage isolates belonged to ST2 which was found to be the prevalent circulating ST in the clinical settings of India and also in global cases. The P. aeruginosa isolates (774 clinical and 06 sewage) segregated into different clades and the study isolates clustered with the clinical isolates from India, East Asia and North America (Fig. 7e). The ST357, one of the dominating STs in the clinical cases of India and East Asia was also observed in sewage. In the case of E. faecium, the isolates (639 clinical and 21 sewage) were quite diverse on the phylogenetic tree (Fig. 7f). The most prevalent ST in sewage isolates was ST80 which was also observed to be prevalent in both Indian and global clinical settings. Although most of the study isolates were positioned adjacent to the clinical isolates from Europe and East Asian countries on the tree, few of the sewage isolates clustered with the Indian clinical isolates. The M. morganii isolates (196 clinical and 33 sewage) were clonal forming distinct clades with most isolates clustering with clinical isolates from East Asia, Europe, Oceania and South America (Supplementary Figs. 2a). The P. mirabilis sewage isolates are segregated to different clades clustering with clinical isolates from North America (Supplementary Fig. 2b).

Genetics of antibiotic resistance traits of MDR bacterial pathogens

The analysis revealed the presence of numerous ARGs encoding resistance to majority classes of antibiotics, which were species-specific as well as common among different bacteria. Morganella, Citrobacter and Aeromonas showed the presence of limited ARGs when compared to the other Gram-negative pathogens such as, Escherichia, Enterobacter, Pseudomonas, Klebsiella and Acinetobacter spp.. The resistance encoded was either through enzymatic inactivation, target alteration or antibiotic efflux (Fig. 8a-c & Supplementary Fig. 3a-c). Majority of ARGs encoded resistance to beta-lactams. The bla_ADC, bla_SHV, bla_OKP, bla_PDC, bla_CPH were specific to A. baumannii, K. pneumoniae, Klebsiella quasipneumoniae, P. aeruginosa, and A. veronii, respectively, with almost all the genomes of these species having the specific ARGs. Although not specific to just M. morganii, bla_DHA was found to be quite prevalent in M. morganii with the presence of one or the other allele in its genome. Except in Enterococcus, the Extended-spectrum beta-lactamase (ESBL) bla_OXA was found in all the species. However, the allele varied from species to species, some being organism-specific, such as bla_OXA−66 and bla_OXA−23 were specific to A. baumannii, bla_OXA−50 to P. aeruginosa, and bla_OXA−12 to A. veronii. Similarly, the carbapenemase bla_NDM−1 was observed in Providencia, Pseudomonas, Klebsiella and Enterobacter, while bla_NDM−5 was observed only in E. coli. The most prevalent beta-lactamase genes in E. coli were the ESBLs bla_TEM (n = 37/69), followed by bla_CTX (n = 23/69). Intriguingly, none of the Enterococcus isolates showed the presence of any beta-lactamase genes (Fig. 8a, Supplementary Table 8d-l).

Fig. 8

Antibiotic resistant genes profile in the sewage isolates

Heatmap depicts the presence (filled) of ARGs in the study isolate, showing resistance to different classes of the antibiotics. The genes could be conferring resistance through a) enzymatic inactivation b. efflux c. target modification.

The aminoglycoside modifying enzymes (AMEs) viz. aph(3’’)-Ib or strA and aph(6)-Id or strB were quite prevalent in Acinetobacter, Pseudomonas, Proteus, Escherichia and Klebsiella. Interestingly, Proteus genomes had more of AMEs compared to other ARGs. While armA was prevalent in A. baumannii, rmts were observed in Pseudomonas and Enterobacterales. The aac(6’)-Ie-aph(2’’)-Ia complex and ant(6)-Ia were observed only in E. faecium. The allele aac(6’)-Ii was specific to the Gram-positive Enterococcus and was present in almost all of the E. faecium (n = 20/21). For resistance to macrolides, the msrE and mphE were prevalent in the non-fermenters such as Acinetobacter, Pseudomonas and Enterobacter spp. The mphA was observed more in E. coli and Klebsiella, while mphB was detected in almost all E. coli (n = 66/69) and msrC was specific to E. faecium and present in almost all of the isolates (n = 19/21). Resistance to tetracyclines was mostly encoded by efflux genes and amongst tet(A) and tet(B) were the most prevalent ones. All the 8 Acinetobacter isolates had the tet(B) and the adeABC and the adeFGH efflux system genes for resistance towards tetracyclines. The tet(E) is specific to Aeromonas, tet(G) to Pseudomonas, tet(J) to Proteus, and tet(L), tetM and tetU to Enterococcus. Several efflux genes, including those encoding multidrug efflux proteins, were identified in the isolates, with many being species-specific. E. coli had the highest number of efflux genes, followed by P. aeruginosa where nearly all isolates carrying Pseudomonas-specific mex, mux, opm, opr, triA, triB and triC genes confer resistance to multiple drug classes (Fig. 8b). Genes providing resistance to trimethoprims (dfr) were harboured by all the species, with dfrA1 being the most prevalent, majorly observed in Providencia, Proteus, Morganella, and Escherichia. Except for Enterococcus, genes encoding resistance to sulfonamides (sul1 and sul2) were observed in all the species. Enterococcus on the other hand, showed resistance to vancomycin (van genes) (Fig. 8c).

Comparative analysis of the genomes from sewage pathogens and those from the global clinical isolates revealed notable similarities (Supplementary Table 7a-j). The ESBL bla_CTX−M−15 was the most prevalent allele of the bla_CTX−M genes in Escherichia, Klebsiella and Morganella in the sewage as well as clinical isolates. While bla_NDM−5 was the most prevalent allele of metallo-β-lactamase in the Escherichia in clinical cases and wastewater, bla_NDM−1 is the most prevalent allele in Klebsiella and Morganella from clinics and wastewater. The cephalosporin resistant bla_DHA−20 was quite prevalent in the clinical and sewage Morganella isolates. However, bla_DHA−1 was almost equally prevalent in wastewater, this was not as frequently observed in the clinical cases. Likewise, the bla_DHA−17 that was the most prevalent beta-lactamase in the Morganella isolates from clinics, this allele was completely absent in the Morganella from wastewater. Strikingly, the prevalence of strA and strB genes in the sewage samples were identified to be 50% more as compared to the clinical samples. While vanA was the prevalent gene responsible for vancomycin resistance in the global isolates followed by vanB gene, the Indian clinical isolates exclusively had only vanA gene. Interestingly, this was also reflected in the sewage isolates and they only carried the vanA gene.

Genetic linkage of antibiotic resistance genes with mobile genetic elements

Many of these ARGs were found to be associated with MGEs such as, transposons, insertion sequences (IS), and miniature inverted-repeat transposable elements (MITEs). A total of 187 out of 305 genomes showed the association of one or more ARGs with MGEs and IS elements were identified to be the most abundant (Fig. 9a). The species with most number of isolates having ARGs linked with either IS elements or transposons are E. coli (n = 61/69), followed by K. pneumoniae (n = 25/33), M. morganii (n = 20/33), P. mirabilis (n = 14/25), A. baumannii (n = 07/08) and P. putida (n = 06/06). The IS6100 was the predominant IS element found in the study isolates. The most common ARGs that were linked with MGEs are aph(3’’)-Ib, aph(6)-Id, ant(3’’)-IIa, the dfr genes and the sul genes (Supplementary Table 8a). While in most of the cases the MGEs linked with these ARGs were different across the species, some ARGs were linked with the same MGE across the species. For example, Tn7 is associated with ant(3’’)-IIa and dfrA1 in P. mirabilis and M. morganii and bla_CTX−M−15 is linked with ISEc9 in E. coli, K. pneumoniae and M. morganii implying inter-species movement of the ARGs by MGEs. In a few isolates of E. coli and M. morganii, the chloramphenicol resistance genes (catI and catB3) were linked with Tn6196. In E. coli, ARGs such as mphA, sul1, bla_DHA−1, qnrB4, aadA5 and the dfr genes were associated with IS6100 but in P. putida, IS6100 is associated with mphE, msrE, aph(3’’)-Ib and aph(6)-Id. In majority of E. coli isolates that showed linkage of strA, strB and bla_TEM−1 with a MGE, the MGE has been IS1133 and/or IS5075 but in A. baumannii, the MGE associated with these two AMEs is ISvsa3. The Gram-positive Enterococcus harboured a completely different set of IS elements and transposons as compared to the Gram-negative species (Fig. 9a).

Fig. 9

Mobile genetic elements associated with antimicrobial resistance genes in the sewage isolates

a) Sankey diagram shows the different ARGs associated with the different MGEs in the different organisms from the sewage. Most ARGs are associated with the IS6100 element in different organisms. E. coli has the highest number of ARGs associated with MGEs. b) Sankey diagram shows the ARGs in different sewage isolates that are linked with plasmid genes. The Inc plasmids are the most prevalent. c) Chord diagram shows the ARGs in different organisms that are found within integrons. The thickness of the lines (connections) indicates the prevalence of that association.

Further, these isolates were screened for the presence of acquired plasmids through the PlasmidFinder database. The database identified the presence of acquired plasmid genes in Klebsiella, Escherichia, Proteus, Enterobacter, Citrobacter and Enterococcus (Supplementary Table 8b). However, very few isolates showed the presence of ARGs in these plasmids. The beta-lactamase genes bla_OXA−181, bla_OXA−232, bla_OXA−484, and bla_CTX−M−15 were associated with the ColKP3 plasmid in E. coli which also harboured many Inc plasmids that carried AMEs such as aph(3’’)-Ib, aph(6)-Id and aph(3’)-Ia. Apart from E. coli, a few isolates of Klebsiella, Citrobacter and Enterococcus showed the presence of ARGs with plasmids. Other ARGs carried on Inc plasmids were sul2, tet(A), bla_CMY−59, and bla_TEM−1 (Fig. 9b). IntegronFinder identified integrons in all the species. In many of the isolates of Acinetobacter, Aeromonas, Citrobacter, Enterobacter, Escherichia, Klebsiella, Morganella, Proteus, Providencia, and Pseudomonas, ARGs were identified within these integrons (Supplementary Table 8c). Most of these ARGs were the AMEs aadA2 and ant(3’’)-IIa, trimethoprim resistant dfr genes and the rifamycin resistant arr-2 gene (Fig. 9c).

Discussion

Sewage is known to be a hotspot for bacterial HGT which serves as a good proxy to survey the diversity of AMR and for identifying the resistant pathogens associated with humans and animals^5–6. Analyzing wastewater is a cutting-edge strategy due to environmental factors like high exposure to surface runoff, direct human or animal contact, or untreated wastewater. Previously, the MDR bacteria confined to frail and debilitated hospital patients, are now spreading in communities posing a growing threat to public health. Analyzing wastewater systems helps in identifying circulating resistant bacteria and track emerging resistance trends, providing an early warning system for potential outbreaks. In this study, we have investigated the AMR spectrum of open drainage systems using both culture-dependent and culture-independent metagenomic approaches.

Antibiotics are widely prescribed and commonly enter the environment through improper drug disposal, excreted waste or municipal dumpsites containing human or animal pharmaceutical waste⁷. Our analysis demonstrated widespread presence of clinically significant antibiotics with tetracycline, amoxicillin and spectinomycin showing particularly high concentrations. The detection frequency and concentration of antibiotics in the wastewater could be influenced by the site, or seasonal variation in the collection of the sample⁸. However, in contrast, antibiotics like ampicillin and polymyxin-B were not detected, possibly due to lower usage or effective removal from the sewage system. This variation in the detection of antibiotics could be attributed to different factors affecting their persistence and distribution, such as their class, molecular weight, and chemical stability with unstable antibiotics degrading faster from the environment. Studies have found various antibiotics in sewage globally^9–12, with ciprofloxacin, azithromycin, and cefalexin prevalent in Europe¹², sulfamethoxazole and trimethoprim common in Africa^13–14, and high levels of fluoroquinolones in Asia^15–16. In India, as per the previous reports wastewater frequently contains sulfamethoxazole, amoxicillin, fluoroquinolones, and trimethoprim ^8,11 which is in accordance with our findings. However, the current study also identified elevated concentrations of different antibiotics which is particularly concerning as they may accelerate the dissemination of ARGs from environment to human.

The dipstick based assay we have standardized in this study demonstrates high specificity and the results were in accordance when compared to the shotgun sequencing data. This nucleic-acid based detection method could be useful for the early detection of ARGs in different settings and the same can also be used for the detection of pathogens directly from the clinical or environmental samples. This assay is rapid and cost effective when compared to shotgun sequencing for epidemiological studies or to monitor the pathogens and associated ARGs in a community.

A large number of taxa, representing a wide phylogenetic span, was uncovered from the sewage samples from community and hospital settings. Interestingly, we discovered notable differences in how seasonality affected certain bacterial taxa. These variations suggest that specific bacterial groups respond differently to seasonal changes, which may influence their abundance, diversity, or pathogen dynamics in different environments^17–19.

Our finding is consistent with previous research that highlights the prominence of beta-lactam resistance in both clinical and environmental contexts²⁰. Research by Li et al., in 2015²¹, highlighted the role of these efflux pumps in both clinical and environmental bacterial isolates, reinforcing efflux mechanisms as a key factor in the resistance landscape. Our identification of 66 efflux genes supports previous studies highlighting efflux pumps as a key mechanism of resistance in bacteria. The varying abundances of ARGs between the community and hospital settings could reflect differences in antibiotic usage patterns and wastewater treatment processes. Kraemer et al. in 2019²² had previously discussed how hospitals typically have higher selective pressures due to intensive antibiotic use, leading to enriched resistance profiles in their effluents. Our findings indicate a significant overlap between sewage and clinical isolates, underscoring the potential role of wastewater environments as reservoirs for drug-resistant bacteria. Similar findings have been reported by others that highlight the presence of clinically relevant strains in environmental sources ^{20, 23}.

The diversity within E. coli isolates, with 22 different sequence types (STs), suggests that sewage harbors significant genetic variation indicating the genetic diversity of MDR E.coli in the study population. Furthermore, the prevalence of K. pneumoniae STs (ST15, ST37, and ST16) in sewage is mirrored by their prevalence in clinical cases across America and Asia indicating a shared lineage that likely facilitates transmission of these strains among clinical and environmental settings. This underscores the critical role of specific STs in the global epidemiology of antibiotic resistance, with these lineages being implicated in major outbreaks worldwide ^24–25.

The widespread presence of ARGs across species highlights the growing threat of XDR pathogens, with prominent challenge of beta-lactam resistance driven by bla_NDM and bla_OXA genes which is in consistence with global data ^20,26. The current study observed differences in ARG profiles among genera, suggesting diverse selective pressure in environment and clinical settings, revealing potential gaps in resistance mechanisms for targeted intervention.

Furthermore, many of these AMR genes were associated with MGEs and the finding that certain ARGs, such as aph(3’’)-Ib and bla_CTX−M−15, are linked with the same MGE across different species implies a significant inter-species transfer of resistance traits, further complicating the landscape of AMR. These findings are consistent with studies demonstrating that MGEs such as insertion sequences, transposons and integrons are critical for the evolution of MDR bacteria ^27–28. However, the specific association of resistance genes with distinct MGEs, such as IS1133 and IS5075 in E. coli or ISvsa3 in A. baumannii, suggests species-specific adaptation strategies that may influence the effectiveness of treatment options.

Our findings indicate that the vanA gene is the dominant driver of vancomycin-resistant Enterococcus (VRE) in India, consistent with global reports ^29–30. In contrast, vanB isolates were less common compared to global data, suggesting geographical or environmental factors influence resistance patterns. Notably, resistance genes like bla_DHA−20 in Morganella and strA and strB for aminoglycoside resistance were more prevalent in sewage isolates than in clinical samples, highlighting a significant environmental reservoir.

In conclusion, this study highlights the pivotal role of sewage as a significant reservoir for MDR and XDR pathogens, facilitating the transfer of ARGs into clinical environments. The findings underscore the influence of antibiotic residues and MGEs in promoting HGT, which fosters bacterial adaptability and resistance. Furthermore, the dipstick-based assay developed in this research provides a rapid, cost-effective, and specific tool for sewage surveillance, allowing for more timely and targeted interventions. The shared MDR lineages and resistant alleles between clinical and environmental settings emphasize the interconnection of these environments and reinforce the necessity of a One Health approach. In light of these insights, the study advocates for the implementation of nationwide sewage monitoring programs and the development of region-specific strategies for outbreak preparedness. Tailoring antibiotic use policies to the unique conditions of each region, rather than relying on generic guidelines, will be essential in mitigating the growing threat of antibiotic resistance.

Methods

Sewage collection and measurement of antibiotic concentration using QTRAP mass spectrometer

Wastewater samples (n = 371) were collected from seven different sites in the city of Faridabad, Haryana (North India) including both hospital (n = 82) and community sewage outlets (n = 250) for a duration of six months from June 2023 to December 2023. Apart from Faridabad, samples (n = 39) were also collected from other parts of India including Assam (Northeast), West Bengal & Jharkhand (Eastern India), Uttar Pradesh, and Uttarakhand (Northern India). About 30 ml of wastewater from each site was collected at least twice a week and was immediately placed on ice and transported to the Functional Genetics Laboratory, THSTI, India on ice for further processing.

The Sciex Exion LC system paired with the QTRAP 6500 + MS was used for LC-MS/MS analysis of antibiotic residues in sewage water. Samples were filtered through a 0.22 µm Millex GV syringe filter (Merck, Darmstadt, Germany) before analysis. Separation was performed on an Acquity UPLC BEH C18 column at 40°C with a gradient elution of 0.1% formic acid in water (A) and acetonitrile (B), ramping B from 15–90% in six minutes at 0.2 mL/min. For electrospray ionization, electrospray voltage used was + 5500 V, with a source temperature of 500°C. Transition ranges were determined by initially injecting each commercially purchased pure antibiotic into the LC-MS/MS. One transition per antibiotic with the sharpest peaks was selected for identification purposes and ANALYST 1.7.2 software was used for data analysis ³¹.

Isolation of sewage metagenomic DNA and screening of antibiotic resistance genes by dipstick based assay

Metagenomic DNA was extracted from sewage samples using the THSTI DNA extraction method³² which includes a series of chemical, physical and mechanical actions to disrupt the cells and release the genomic DNA. The precipitated genomic DNA was washed and resuspended in nuclease-free water followed by resolving the sample in 0.8% agarose and purity check.

The purified metagenomic DNA was used as a template to detect the presence and prevalence of ARGs in sewage water. A total of sixteen emerging resistant determinants (Supplementary Table 2) coding resistance to different drug classes routinely used in hospital settings or are highly prevalent among bacterial pathogens were selected. The designed tagged and biotinylated set of primers specific for each resistant allele was utilized to standardize the multiplex PCR assays. Further, the m-PCR assays were translated into C-PAS (4) dipstick test strips (TBA Co., Ltd) using the amplicons along with the developing solutions mixing with the streptavidin-coated blue-coloured latex beads. The immersion of the dipstick in the resulting solution generates a blue colour line indicating the detection of target genes.

16S rRNA gene and shotgun sequencing of sewage metagenome

The metagenomic DNA of the microbiome was amplified for the targeted metagenomics of the 16S rRNA specific V3 to V4 region. The nucleoside nomenclature of the overhang (shown in italics) holding primers for targeting the amplification sites are: Forward Primer: 5'TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACGGGNGGCWGCAG3', Reverse Primer: 5'GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTACHVGGGTATCTAATCC3'. The amplicons were cleaned up with the Agencourt AMPure XP paramagnetic beads (Beckman Coulter). Samples were subjected for the library preparation with the Nextera XT Index Kit V2 (Illumina). The cleaned-up libraries were quantified with Qubit Flex Fluorometer (invitrogen) and verified for quality analysis with TapeStation 4200 (Agilent) automated electrophoresis system. Qualified libraries have been normalized to 4nM concentration and a final concentration of 700pM pooled library was sequenced in the Nextseq 2000 sequencing platform (Illumina) with paired end 2X300 cycle chemistry.

For the shotgun sequencing the metagenomic DNA samples were quantified with the Qubit dsDNA BR Assay Kits and HS Assay kits (invitrogen) in the Qubit Flex Fluorometer. A good quantity of samples with 300ng were tagmented to short fragments using the bead-linked transposome (Illumina). The library was prepared using the washed fragments with the Illumina DNA prep kit and IDT-Illumina DNA/RNA UD indexes (Illumina). The library was cleaned up following left-side and right-side size selection double clean up method with the Agencourt AMPure XP paramagnetic beads (Beckman Coulter) and the library quantity was determined as described earlier. The TapeStation 4200 (Agilent) has been used to check the bell-shaped curve formation in the electropherograms of the libraries. Presence of any impurities was cleaned up using the SPRIselect beads (Beckman Coulter). The libraries were normalized to 4nM concentration and pooled for sequencing at a final concentration of 750pM in the Nextseq2000 sequencing platform (Illumina). The sequencing chemistry used the paired end 2X150 cycles chemistry in the P2 flow cell (Illumina) and the data generated has been submitted to NCBI (16s metagenome: Bioproject ID: PRJNA1130090; Shotgun data: PRJNA1135695) and the IBDC (Accession number: INRP000179) database.

Analysis of bacterial diversity and antibiotic resistance genes in the sewage metagenome

Demultiplexed FASTQ files for Read 1 (R1.fastq) and Read 2 (R2.fastq) of each sample were subjected to initial quality control based on the Fastp (version 0.20.1) reports generated. The trimmed reads were processed using DADA2 package³³ (Version 1.16) in R (Version 4.1.1) to generate amplicon sequence variant (ASV) generation. For species-level classification, the NCBI 16S Microbial database was used with blaSTn for alignment³⁴. Only representative sequences with > 98% sequence identity in BLAST were selected for species-level annotation. For shotgun data, raw reads were processed and trimmed using as described above and the resulting reads were subsequently aligned to the human reference database using HISAT2. The reads that remained unaligned after HISAT2 were concurrently subjected to taxonomic classification using Kraken2³⁵ and assembled with SPAdes (version 3.13.0) in meta mode³⁶. Bacterial reads were extracted selectively via a custom bash script after Kraken2. The bacterial diversity within the sewage metagenome was assessed using the phyloseq package in R. Metagenome binning for each sample was carried out with MetaBAT2³⁷. Metagenome-assembled genomes (MAGs) were screened for ARGs using the abricate tool with the CARD database, applying default settings (> 80% identity, > 80% coverage). The contigs with CARD hits were taxonomically classified with Kraken 2, and taxa were evaluated from species to phylum levels. The contigs/scaffolds that met the following criteria were selected for network analysis of bacteria and ARGs: (1) classified as belonging to the superkingdom Bacteria by Kraken 2, (2) contained an ARG as identified by the CARD tool, and (3) had a Kraken 2 taxonomic classification at the species level. Bacteria and ARG nodes connected in the network were visualized using Cytoscape. Additionally, the MAGs were queried against KEGG Orthology (KO) IDs, Clusters of Orthologous Groups (COG), and Gene Ontology (GO) using EggNOG-mapper, with parameters set for > 90% identity and > 80% coverage³⁸.

Isolation, identification and cultivation of drug-resistant aerobic and anaerobic bacterial species

The sample resuspended in enriched medium, was inoculated into Mueller-Hinton and Brain Heart Infusion (BHI) broth, under selective pressure from five different antibiotics: ampicillin (100µg/ml), streptomycin (100µg/ml), nalidixic acid (10µg/ml), azithromycin (15µg/ml), and tetracycline (10µg/ml). Samples demonstrating optimal growth were subsequently processed for culturomics analysis. This involved further culturing on MacConkey and BHI agar to isolate both aerobic and anaerobic organisms. Isolates were identified using 16S rDNA Sanger sequencing, and the results were analyzed in the NCBI database.

Antimicrobial susceptibility testing of bacterial isolates

The susceptibility profiles of a total of 681 aerobic and facultative anaerobic isolates were assessed using the Kirby-Bauer disc diffusion method against a panel of 32 commercially available antibiotic discs. The interpretation of results adhered to the latest guidelines from the Clinical and Laboratory Standards Institute (CLSI) and The European Committee on Antimicrobial Susceptibility Testing (EUCAST)^39–40. Reference strains, including E. coli ATCC 25922, A. baumannii ATCC 17978, and K. pneumoniae ATCC 700603, were employed to ensure the accuracy of the experiments. The specific details on antibiotic concentrations and bacterial isolates tested are provided in the supplementary Table 5.

Whole genome sequencing and bioinformatics analysis

The bacterial isolates confirmed with the Sanger sequencing (Applied Biosystems Genetic Analyzer 3500) were subjected for the whole genome sequencing (WGS) followed by tagmentation with in vitro transposition. The transposed fragments readily bind to the index adapters, Nextera XT Index kit V2 from Illumina and the Nextera XT library preparation kit. After indexing, the library pool was sequenced at 750pM following the methods as described for shotgun sequencing. The raw files were processed using the seqtk tool v1.3-r106 and then merged using PEAR (Paired-End read merger) software. The SPAdes pipeline was employed to assemble the genome from the cleaned paired-end reads³⁶. The quality of the assembled genomes was evaluated using CheckM (version v1.1.3)⁴¹. Only genomes meeting the criteria of > 90% completeness and < 5% contamination were included for subsystem analysis. Annotation was carried out using PROKKA v1.14.6. The whole genome sequences of the bacterial isolates are submitted to the NCBI (Bioproject ID: PRJNA1128678) and the IBDC (Accession number: INRP000179) database.

Multilocus sequence typing

The sequence type (ST) of the genomes was identified using the MLST finder (https://github.com/tseemann/mlst). The Achtman scheme was used for E. coli (adk, fumC, gyrB, icd, mdh, purA, recA), the Pasteur scheme was used for A. baumannii (cpn60, fusA, gltA, pyrG, recA, rplB, rpoB) and the PubMLST scheme was used for K. pneumoniae (phoE, mdh, tonB, gapA, rpoB, pgi, infB), P. aeruginosa (acsA, areE, guaA, mutL, nuoD, ppsA, trpE) and E. faecium (adk, atpA, ddl, gdh, gyd, pstS, purK)⁴².

Phylogenetic analysis

The genomes of the study isolates were mapped against the consensus sequences of the multiple copies of the 5S, 16S, 23S rRNA genes of the NCBI reference sequences using Snippy v4.6.0 (https://github.com/tseemann/snippy). For global phylogeographical analysis, the whole genome assemblies of A. baumannii (n = 689), E. coli (n = 2700), K. pneumoniae (n = 1281), P. aeruginosa (n = 774), M. morganii (n = 196), P. mirabilis (n = 513), P. rettgeri (n = 132), C. werkmanii (n = 22), and E. faecium (n = 639) from clinical specimen (Homo sapiens) submitted during 2019–2024 from across the globe were downloaded from NCBI. Only the genomes that passed CheckM criteria with > 90% completeness and < 5% contamination were included for the analysis⁴¹. Alignment of core gene SNPs was generated using Snippy v4.6.0. The maximum likelihood phylogenetic trees were constructed by RAxML v8.2.12 using the GTRGAMMA model⁴³. The trees were visualized using iTOL and annotated with metadata such as year of isolation, geographic location, and sample source⁴⁴. Clermontyper v20.03 was used to type the E. coli genomes into phylogroups⁴⁵.

Identification and analysis of antibiotic resistance genes and mobile genetic elements

The CARD database and the PlasmidFinder in the ABRicate tool v1.0.1 were used to predict the antimicrobial resistance genes and plasmid genes, respectively, in the isolates (https://github.com/tseemann/abricate). For all the databases, a filtering criterion of 90% identity and 80% gene coverage was used. The MGEs were identified using the MobileElementFinder v1.1.2 that integrates data from ICEberg, ISFinder and the Tn repository⁴⁶. The Integron Finder v2.0.5 was used to identify the presence of integrons in the isolates⁴⁷.

Statistical analysis

Alpha diversity indices were computed using the ‘estimate_richness’ function from the phyloseq R package. Beta diversity analyses were carried out with the vegan and phyloseq R packages. Differences in alpha diversity between groups were assessed using the non-parametric Wilcoxon rank-sum test (two-tailed) with the ‘wilcox.test’ function, and boxplots were generated using the ggpubr R package. Statistical analysis and microbial composition comparisons among groups were performed using the ampvis2 and ggplot2 R packages. Core taxa were identified with the eulerr, microbiome, and microbiome utilities R packages. A p-value of 0.05 or less was considered statistically significant.

Data availability

All genome sequences data of 16s metagenomics, Shotgun sequencing and WGS data generated are deposited in the National Center for Biotechnology Information (https://www.ncbi.nlm.nih.gov/), and Indian Biological Data Centre (https://ibdc.dbtindia.gov.in). The submission IDs are included in the specific Materials and Methods section. Also, all the genome sequence data used in the study are available in the NCBI database.

Acknowledgements

We extend our sincere gratitude to Prof. Pramod K Garg and Prof. Ganesan Karthikeyan for their generous support to complete this study. Authors also acknowledged Prof. Sudhansu Vrati, Prof. Arindam Maitra, and Dr. Onkar Nath Tiwari for their support to the execution of this research. This study was funded by the Department of Biotechnology, Government of India, under the surveillance project titled “Genomic Surveillance for SARS-CoV-2 in India: Indian SARS-CoV-2 Genomics Consortium (INSACOG) - Phase II: Sewage Surveillance” (Ref. RAD/22017/19/2022-KGD-DBT, dated 29/12/2022).

Author contributions

B.D. conceptualized and designed the research. D.P, D.T, R.S.K, V.G, L.N, P.J, P.K, J.S, S.K, C.B, K.K, S.B, S.L., S.T, R.K, P.B, M.B., and Y.K performed experiments. S.M, and S.K.B, contributed to the investigation and sample collection. B.D., D.P and A.M contributed to the development of dipstick based detection of ARGs. D.P, D.T, K.R.S, N.W, and B.D contributed to formal analysis, validation and writing the original draft. All the authors approved the final version of the manuscript.

Competing interests:

The authors have declared that no competing interests exist.

Additional information

Supplementary information

The supplementary materials for this study is available at Supplementary Table S1 to S8.

Correspondence and request for materials

should be addressed to Dr Bhabatosh Das, Professor, THSTI.

Legends

Supplementary Fig. 1: Diversity measure and xenobiotic pathways in the community and hospital setting of Faridabad

a) Diversity measures i) Shannon index (Wilcoxon, p-value = 0.0041) ii) Chao1 (Wilcoxon, p-value = 0.0067) iii) Simpson (Wilcoxon, p-value = 0.0024). The box indicates the interquartile range (IQR). The median value is represented as a line within the box and whiskers extend to the extreme value i.e. within 1.5 × IQR. b) Bar plot showing the list of pathways likely involved in the degradation of xenobiotic pollutants. 17 KEGG categories associated with various xenobiotic pollutants are listed with Benzoate being the most prevalent pollutant in both the sites.

Supplementary Fig. 2: Phylogenomic trees based on core-genome SNPs showing the genetic relatedness of sewage isolates with clinical samples

a. M. morganii b. P. mirabilis c. P. rettgeri d. C. werkmanii.

Supplementary Fig. 3: Heatmap shows the percent genomes of the studied genus having ARGs against different antibiotic classes

. The resistance could be through a. enzymatic inactivation b. efflux c. target modification.

Supplementary Information

Table S1

Concentration distribution of various antibiotic residues in wastewater samples. The various sheets in the Excel file provide data on the concentrations of different antibiotics detected in sewage water viz. a. Ampicillin; b. Azithromycin; c. Nalidixic acid; d. Amoxicillin; e. Kanamycin; f. Erythromycin; g. Rifampicin; h. Tetracycline; i. Spectinomycin; j. Polymyxin-B; k. Colistin.

Table S2

Detection of antibiotic resistant genes using dipstick based assay. Four sets of multiplex PCR assay were performed to target 16 ARGs, in Group 1 (bla_NDM, bla_CTX−M, bla_KPC, bla_TEM); Group 2(cat, aac, sul2, oqxA); Group 3 (bla_OXA, tet, aadA, sul1); Group 4 (ermB, dfrA, qnrS, mphA) were targeted and detected using dipstick assay.

Table S3

Read information from 16S sequencing samples processed by DADA2 pipeline

Table S4

Distribution of sewage resistome data for a) ARGs in community and hospital settings, b) efflux genes in community and hospital settings

Table S5

Susceptibility profile of the sewage isolates a) Isolates were checked towards different classes of antibiotics using Kirby-Bauer disc diffusion method and were categorized as resistant, susceptible, and intermediate. The diameter of each of the antibiotics are also mentioned along with the antibiotics which are intrinsically resistant (IR) to the isolate. b) Screening of vancomycin-resistant Enterococcus faecium isolates (VRE) among the sewage isolates.

Table S6

Percentage of resistance of the isolates a) Overall count of resistance towards antibiotics b) Total percentage of the organism resistant towards different classes of antibiotics

Table S7

Sample information of the sewage whole genome sequences and the NCBI-downloaded clinical sequences. a) Information of whole genome sequences from sewage isolates. The sheets b-k are the metadata of the clinical isolates downloaded from the NCBI b) E. coli c) K. pneumoniae d) A. baumannii e) P. aeruginosa f) E. faecium g) P. mirabilis h) M. morganii i) P. rettgeri k) C. werkmanii

Table S8

Genomic analysis information on the ARGs, MGEs, plasmids and integrons. a) MGEs associated with ARGs in sewage isolates b) ARGs associated with plasmid genes in sewage isolates c) Integrons linked with ARGs in sewage isolates. The sheets d-l represents the presence of ARGs in different clinical isolates d) A. baumannii e) E. coli f) K. pneumoniae g) P. aeruginosa h) M. morganii i) P. mirabilis j) P. rettgeri k) C. werkmanii l) E. faecium.

Electronic Supplementary Material

Below is the link to the electronic supplementary material

Supplementary Material 1

Supplementary Material 2

Supplementary Material 3

References

Naghavi M et al (2024) Global burden of bacterial antimicrobial resistance 1990–2021: a systematic analysis with forecasts to 2050. Lancet 404:1199–1226

Micoli F, Bagnoli F, Rappuoli R, Serruto D (2021) The role of vaccines in combating antimicrobial resistance. Nat Rev Microbiol 19:287–302

Bottery MJ, Pitchford JW, Friman VP (2021) Ecology and evolution of antimicrobial resistance in bacterial communities. ISME J 15:939–948

FAO-OIE-WHO One Health approach to AMR mitigation and safer food in the Asia-Pacific-Region, Regional Office for South-East Asia; License: CC BY-NC-SA 3.0 IGO (2022)

Munk P et al (2022) Genomic analysis of sewage from 101 countries reveals global landscape of antimicrobial resistance. Nat Commun 13:7251

Sambaza SS, Naicker N (2023) Contribution of wastewater to antimicrobial resistance: A review article. J Glob Antimicrob Resist. ; 34: 23–29 (2023)

Rayan RR (2023) Pharmaceutical effluent evokes superbugs in the environment: A call to action. Biosaf Health 5:363–371

Rathinavelu S et al (2024) Mapping the scarcity of data on antibiotics in natural and engineered water environments across India. Front Antibiot 3:1337261

Hu J, Zhou J, Zhou S, Wu P, Tsang YF (2018) Occurrence and fate of antibiotics in a wastewater treatment plant and their biological effects on receiving waters in Guizhou. Process Saf Environ Prot 113:483–490

10.

Zainab SM, Junaid M, Xu N, Malik RN (2020) Antibiotics and antibiotic resistant genes (ARGs) in groundwater: a global review on dissemination, sources, interactions, environ- mental and human health risks. Water Res 187:116455

11.

Arun S et al (2022) Antibiotics in sewage treatment plants, receiving water bodies and groundwater of Chennai city and the suburb, South India: Occurrence, removal efficiencies, and risk assessment. Sci Total Environ 10:851:158195

12.

Sara Rodriguez-Mozaz S et al (2020) Antibiotic residues in final effluents of European wastewater treatment plants and their impact on the aquatic environment. Environ Int 140:105733

13.

Ncube S, Nuapia YB, Chimuka L, Madikizela LM, Etale A (2021) Trace Detection and Quantitation of Antibiotics in a South African Stream Receiving Wastewater Effluents and Municipal Dumpsite Leachates. Front Environ Sci 9:733065

14.

Madikizela LM, Ncube S, Chimuka L, Analysis (2020) Occurrence and Removal of Pharmaceuticals in African Water Resources: A Current Status. J Environ Manage 253:109741

15.

Shi Y et al (2020) Antibiotics in Wastewater from Multiple Sources and Surface Water of the Yangtze River in Chongqing in China. Environ Monit Assess. 192

16.

Zafar R, Bashir S, Nabi D, Arshad M (2021) Occurrence and Quantication of Prevalent Antibiotics in Wastewater Samples from Rawalpindi and Islamabad, Pakistan. Sci Total Environ 764:142596

17.

Kang XH, Leng Y, Zeng OMM, X.Y., Li SW (2020) The seasonal changes of core bacterial community decide sewage purification in sub-plateau municipal sewage treatment plants. Bioprocess Biosyst Eng 43:1609–1617

18.

He Y et al (2021) Microbial Community Structure of Waste Water Treatment Plants in Different Seasons. Huan Jing Ke Xue 8:42: 1488–1495

19.

Hou L, Zhou Q, Wu Q, Gu Q, Sun M, Zhang J (2018) Spatiotemporal changes in bacterial community and microbial activity in a full-scale drinking water treatment plant. Sci Total Environ 1:625: 449–459

20.

De Angelis G, Del Giacomo P, Posteraro B, Sanguinetti M, Tumbarello M (2020) Molecular Mechanisms, Epidemiology, and Clinical Importance of β-Lactam Resistance in Enterobacteriaceae. Int J Mol Sci 18:21:5090

21.

Li XZ, Plesiat P, Nikaido H (2015) The challenge of efflux-mediated antibiotic resistance in Gram-negative bacteria. Clin Microbiol Rev 28:337–418

22.

Kraemer SA, Ramachandran A, Perron GG (2019) Antibiotic Pollution in the Environment: From Microbial Ecology to Public Policy. Microorganisms 22:7:180

23.

Fouz N et al (2020) The Contribution of Wastewater to the Transmission of Antimicrobial Resistance in the Environment: Implications of Mass Gathering Settings. Trop Med Infect Dis 5:33

24.

Gorrie CL et al (2022) Genomic dissection of Klebsiella pneumoniae infections in hospital patients reveals insights into an opportunistic pathogen. Nat Commun 13:3017

25.

Zhu J et al (2016) Outbreak of NDM-1-producing Klebsiella pneumoniae ST76 and ST37 isolates in neonates. Eur J Clin Microbiol Infect Dis 35:611–618

26.

Nordmann P, Naas T, Poirel L (2011) Global spread of carbapenemase-producing Enterobacteriaceae. Emerg Infect Dis 17:1791–1798

27.

Jeon JH, Jang KM, Lee JH, Kang LW, Lee SH (2023) Transmission of antibiotic resistance genes through mobile genetic elements in Acinetobacter baumannii and gene-transfer prevention. Sci Total Environ 857:159497

28.

Razavi M, Kristiansson E, Flach CF, Larsson (2020) D.G.J. The association between Insertion Sequences and Antibiotic Resistance Genes. mSphere, 5:e00418-20

29.

Azzam A et al (2023) Prevalence of vancomycin-resistant enterococci (VRE) in Egypt (2010–2022): a systematic review and meta-analysis. J Egypt Public Health Assoc 98:8

30.

Werner G, Klare I, Fleige C, Witte W (2008) Increasing rates of vancomycin resistance among Enterococcus faecium isolated from German hospitals between 2004 and 2006 are due to wide clonal dissemination of vancomycin-resistant enterococci and horizontal spread of vanA clusters. Int J Med Microbiol 298:515–527

31.

Chan CL et al (2022) A Universal LC-MS/MS Method for Simultaneous Detection of Antibiotic Residues in Animal and Environmental Samples. Antibiot (Basel) 24:11: 845

32.

Bag S et al (2016) An Improved Method for High Quality Metagenomics DNA Extraction from Human and Environmental Samples. Sci Rep 6:26775

33.

Callahan BJ et al (2016) DADA2: High-resolution sample inference from Illumina amplicon data. Nat Methods 13:581

34.

Johnson M et al (2008) NCBI blaST: a better web interface. Nucleic Acids Res 36:W5

35.

Wood DE, Lu J, Langmead B (2019) Improved metagenomic analysis with Kraken 2. Genome Biol 20:257

36.

Prjibelski A, Antipov D, Meleshko D, Lapidus A, Korobeynikov A (2020) Using SPAdes De Novo Assembler. Curr Protoc Bioinform 70:e102

37.

Kang DD et al (2019) MetaBAT 2: an adaptive binning algorithm for robust and efficient genome reconstruction from metagenome assemblies. Peer J 26:7: e7359

38.

Cantalapiedra CP, Hernandez-Plaza A, Letunic I, Bork P, Huerta-Cepas (2021) J. eggNOG-mapper v2: Functional Annotation, Orthology Assignments, and Domain Prediction at the Metagenomic Scale. Mol Biol Evol 38:5825–5829

39.

CLSI (2021) Performance Standards for Antimicrobial Susceptibility Testing. 31st ed. CLSI supplement M100. Clinical and Laboratory Standards Institute

40.

The European Committee on Antimicrobial Susceptibility Testing (2021) Breakpoint tables for interpretation of MICs and zone diameters. Version 11.0

41.

Parks DH, Imelfort M, Skennerton CT, Hugenholtz P, Tyson GW (2015) CheckM: Assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res 25:1043–1055

42.

Jolley KA, Bray JE, Maiden MCJ (2018) Open-access bacterial population genomics: BIGSdb software, the PubMLST.org website and their applications. Wellcome Open Res 24:3: 124

43.

Stamatakis A (2014) RAxML version 8: A tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30:1312–1313

44.

Letunic I, Bork P (2007) Interactive Tree Of Life (iTOL): An online tool for phylogenetic tree display and annotation. Bioinformatics 23:127–138

45.

Beghain J, Bridier-Nahmias A, Le Nagard H, Denamur E, Clermont O (2018) ClermonTyping: an easy-to-use and accurate in silico method for Escherichia genus strain phylotyping. Microb Genom 4(7):e000192

46.

Neron et al (2022) IntegronFinder 2.0: Identification and Analysis of Integrons across Bacteria, with a Focus on Antibiotic Resistance in Klebsiella. Microorganisms 10:4: 700

47.

Johansson MHK et al (2021) Detection of mobile genetic elements associated with antibiotic resistance in Salmonella enterica using a newly developed web tool: Mobile Element Finder. J Antimicrob Chemother 76:101–109

Yes

Abstract

The emergence and spread of antibiotic resistance (AMR) in clinically important bacterial pathogens severely compromised the effectiveness of commonly used antibiotics in healthcare. Acquisition and transmission of AMR genes (ARGs) are often facilitated by sublethal concentrations of antibiotics in a microbially dense environments. In this study, we used sewage samples (n=371) from six Indian cities between June and December 2023 to assess the concentration of eleven antibiotics, microbial diversity, and ARG richness. Our findings revealed the presence of antibiotics from seven drug classes and over 2000 bacterial ASVs. Metagenomic (n=220) and isolated genome sequences (n=305) of aerobic and anaerobic bacterial species identified 82 ARGs associated with 80 mobile genetic elements (MGEs). These MGEs were predominantly found in multidrug-resistant bacterial pathogens. Comparative core genome analysis (n=7166) of multi drug-resistant (MDR) bacterial isolates showed a strong genetic similarity between sewage-derived strains and clinical pathogens. Our results highlight sewage as a significant reservoir for ARGs, where genetic exchanges occur, facilitating the evolution and spread of AMR pathogens within both community and healthcare settings. Additionally, the dipstick-based assay developed for detection of clinically important ARGs in the present study could be employed for sewage surveillance in low resource settings for better understanding of resistance prevalence.