Urban monitoring of antimicrobial resistance during a COVID-19 surge through wastewater surveillance

During the early phase of the COVID-19 pandemic, infected patients presented with symptoms similar to bacterial pneumonias and were treated with antibiotics before confirmation of a bacterial or fungal co-infection. We reasoned that wastewater surveillance could reveal potential relationships between reduced antimicrobial stewardship, specifically misprescribing antibiotics to treat viral infections, and the occurrence of antimicrobial resistance (AMR) in an urban community. Here, we analyzed microbial communities and AMR profiles in sewage samples from a wastewater treatment plant (WWTP) and a community shelter in Las Vegas, Nevada during a COVID-19 surge in December 2020. Using a respiratory pathogen and AMR enrichment next-generation sequencing panel, we identified four major phyla in the wastewater, including Actinobacteria, Firmicutes, Bacteroidetes and Proteobacteria. Consistent with antibiotics that were reportedly used to treat COVID-19 infections (e.g., fluoroquinolones and beta-lactams), we also measured a significant spike in corresponding AMR genes in the wastewater samples. AMR genes associated with colistin resistance (mcr genes) were also identified exclusively at the WWTP, suggesting that multidrug resistant bacterial infections were being treated during this time. We next compared the Las Vegas sewage data to local 2018–2019 antibiograms, which are antimicrobial susceptibility profile reports about common clinical pathogens. Similar to the discovery of higher levels of beta-lactamase resistance genes in sewage during 2020, beta-lactam antibiotics accounted for 51 ± 3 % of reported antibiotics used in antimicrobial susceptibility tests of 2018–2019 clinical isolates. Our data highlight how wastewater-based epidemiology (WBE) can be leveraged to complement more traditional surveillance efforts by providing community-level data to help identify current and emerging AMR threats.


Introduction
COVID-19, which is caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), resulted in a global pandemic that overwhelmed healthcare systems (Hu et al., 2021;Liu et al., 2020). Early in the pandemic, infected patients exhibited symptoms similar to bacterial pneumonias and were erroneously treated with antibiotics that are ineffective against viral respiratory infections (Nori et al., 2021;Rawson et al., 2020;Lardaro et al., 2021). This adversely impacted antimicrobial stewardship programs that were designed to reduce consequences from excess and inappropriate use. While the occurrence and spread of AMR can be influenced by multiple factors, it is important to characterize the role of the COVID-19 pandemic in potentially compromising these antimicrobial stewardship programs (CDC, 2022).
With this long-term goal in mind, WBE has been utilized to monitor community AMR patterns (Hendriksen et al., 2019;Su et al., 2017;Aarestrup and Woolhouse, 2020;Majeed et al., 2021;Pazda et al., 2019;Charles et al., 2022). Current methods of AMR WBE include the use of quantitative polymerase chain reaction (qPCR) assays to quantify a subset of AMR genes and the use of metagenomic sequencing to characterize both the population of microorganisms and the AMR genes they harbor (Gaviria-Figueroa et al., 2019;Alexander et al., 2020;Guo et al., 2017;Raza et al., 2021;Pärnänen et al., 2019). The advantage of qPCR detection over metagenomics sequencing is often reduced cost, while a disadvantage is the lower throughput of genes detected. Metagenomic sequencing overcomes the gene discovery limitation by identifying an assortment of AMR genes present in a sample, commonly called a resistome, but this typically requires large sequencing depth or target enrichment sequencing to reduce or eliminate non-target reads (Lanza et al., 2018;Guitor et al., 2019;Li et al., 2022). In order to characterize AMR resistomes with low output sequencers, target enrichment sequencing is often necessary and can be accomplished using commercially available target enrichment next-generation sequencing (NGS) panels designed for clinical samples. In clinical settings, the results from cultured-based AMR assays are periodically summarized into antibiograms, which are empirical reports of pathogenic bacterial susceptibility/resistance to a panel of antibiotics. Specifically, bacterial isolates from infected patients are exposed to relevant antibiotics at standard minimum inhibitory concentrations (MICs) to characterize their level of susceptibility/resistance. Healthcare providers then share and use this information on a local scale to provide empiric antimicrobial therapy, which allows clinicians to more strategically prescribe and administer antibiotics without patient-specific culture results, and on a larger scale to combat AMR through antibiotic stewardship (Abbas and Stevens, 2018). Traditional AMR surveillance depends on these time-consuming, culture-based methods that are compiled one patient at a time. However, coupling 'pooled' samples from WBE with emerging NGS technologies offers a potentially more rapid and robust complement for AMR surveillance at the community level.
Here, we investigated the wastewater profile of microorganisms and AMR markers during a COVID-19 surge from November 2020 through January 2021. Sewage samples were collected from a wastewater treatment plant (WWTP) (i.e., community-level surveillance) and from a homeless/rehabilitation shelter (i.e., facility-level surveillance) in the Las Vegas metropolitan area (Nevada, USA). Microbial community profiling was accomplished utilizing a novel respiratory pathogen and AMR enrichment NGS panel that was originally designed to examine clinical samples. The longitudinal analysis of wastewater samples revealed the value of WBE in tracking the evolution of AMR genes, particularly when healthcare resources are strained. This study demonstrates that NGS can be used to complement culture-based data to provide critical insights into AMR in urban environments.

Sample collection, extraction of nucleic acids, and quantification of SARS-CoV-2
Wastewater samples were collected weekly using an autosampler located at the headworks of a WWTP and directly from a manhole that was selected to isolate sewage contributions from a homeless/rehabilitation shelter. The WWTP serves approximately 86,000 individuals and treats an average daily flow of 19 megaliters per day (MLD). Samples collected at the WWTP were 24-h flow-weighted composites of raw influent, and samples were collected at~10:00 am every Monday morning between November 2, 2020 and January 25, 2021. For the shelter manhole sample, flow was not determined at the time of collection, but the shelter was serving fewer than 100 individuals on the collection days (the exact population was not provided by the shelter for privacy reasons). Grab sampling from the shelter manhole was conducted every Monday morning at~6:00 am during a four-week span from November 23, 2020 through December 14, 2020. Sampling was conducted with an autoclaved bottle connected to a Nasco swing sampler and involved collection of three 150-mL grab samples spaced 5 min apart, which were subsequently combined to generate a manual composite (Vo et al., 2022a;Gerrity et al., 2021).

Library preparation and sequencing
Sequencing libraries were constructed using the Respiratory Pathogen ID/AMR Enrichment Panel Kit (RPIP) (Illumina, San Diego, CA, USA) per manufacturer's instructions. The enrichment panel targets >280 respiratory pathogens, including SARS-CoV-2 and other viruses, bacteria, and fungi, in addition to >1200 AMR markers. More than 10 ng of co-extracted RNA/ DNA (i.e., from the PureLink Viral RNA/DNA Mini Kit) was processed for cDNA synthesis using reagents from the RPIP library kit. Libraries were sequenced using an Illumina NextSeq 500 sequencer with a high-output v2.5 (300 cycles) flow cell (Vo et al., 2022a;Vo et al., 2022b). Raw fastq files are available on the National Center for Biotechnology Information website under Bioproject PRJNA833122.

Analysis of local community antibiograms
Publicly available outpatient antibiograms for Clark County, Nevada were obtained for the years 2018-2019 from the Southern Nevada Health District (SNHD) (https://www.southernnevadahealthdistrict.org/clarkcounty-antibiogram/, 2021). Reported antibiotic susceptibilities were grouped according to their antibiotic class and averaged to give a general susceptibility percentage for that class.

Statistical analyses
The Bracken output was used to rarefy the data, and rows with a count <100 were filtered out (Fig. S1). The remaining dataset was used with the vegdist function in R using the Bray method to calculate dissimilarity indices. The Chao1 index was used to measure alpha diversity, while the Bray indices were used for Beta-diversity analysis through non-metric multidimensional scaling. Kruskal-Wallis rank sum tests on abundance of AMR genes belonging to either beta-lactam, fluoroquinolone, or macrolide/ lincosamide/streptogramin (MLS) class were performed to test for statistically significant differences for each sample site. Pairwise Wilcoxon rank sum tests using Benjamini-Hochberg setting for p-value adjustment were performed to test for statistically significant differences in AMR gene class abundance based on month for each sampling site. Pairwise comparisons with a p-value < 0.05 were reported as significant.

Human subjects statement
The University of Nevada Las Vegas Institutional Review Board (IRB) reviewed this project and determined it to be exempt from human subject research according to federal regulations and university policy.

Identification of bacteria in sewage and taxonomic analysis
We first evaluated the relative microbial composition of the sewage samples during November 2020 to January 2021. Sewage samples showed the prevalence of four major phyla, specifically Actinobacteria, Firmicutes, Bacteroidetes and Proteobacteria (Fig. 1A), while dominant taxonomic orders varied based on site and month of sampling (Fig. 1B). For example, Bacteroidales and Campylobacterales were the most abundant orders for the WWTP, while Bacteroidales and Eubacteriales were the most abundant orders for the shelter. Alpha diversity metrics indicated subtle variation in diversity between samples. For example, the highest diversity was observed in a sample from the WWTP, while the lowest diversity was observed in a sample from the shelter; this was also consistent in terms of phylum richness, which was greater for the WWTP than the facility-specific samples (Fig. 1C). These observations could be due to the more complex composition of the community-level wastewater sample, which also traveled longer distances through the sewer collection system. We next examined the non-metric multidimensional scaling (NMDS) of beta diversity (Bray-Curtis), which allows us to measure dissimilarity across samples. Consistent with expectations, our data indicate microbial clustering based on sample location (WWTP vs. shelter) but also greater dissimilarity across the shelter samples (Fig. 1D). These results indicate that microbial community structure in sewage collected at the facility level is significantly different than sewage collected at the community level (i.e., WWTP). Specifically, facility-level sewage may be more likely to exhibit temporal variability and be dominated by fewer taxa. Differences in sample type (i.e., composite vs. grab) may have also contributed to these differences in microbial community structure, particularly for the grab samples at the facility level.

AMR gene detection and abundance relationship with COVID-19 cases
The analysis of sequencing data from the WWTP revealed that AMR gene enrichment accounted for 1.2 % to 12.7 % of the total reads per sample (Table S2), and the abundances and types of AMR classes varied depending on the sampling date and location (Fig. 2). Tetracycline resistance was the most abundant AMR gene class in the WWTP (63 ± 40 %), followed by MLS (28 ± 8 %) (Fig. S2A). Additionally, out of 238 AMR genes detected, the following 10 AMR genes were the most abundant in WWTP samples from November 2021 to January 2022 (listed in decreasing order): tetQ, ermF, tetX, msrE, tetW, mphD, ermG, tetP, mefA, and ermB.
Analysis of the shelter manhole samples in November/December revealed monotonically increasing SARS-CoV-2 levels in wastewater (Fig. 3D), thereby suggesting a COVID-19 outbreak in the facility. A similar wastewater/case surge was observed in the target WWTP for this study and also the WWTP serving the shelter facility (159 MLD and 757,000 people in sewershed; not monitored for AMR in the current study) (Fig. S3). Upon analyzing AMR classes, we identified the most abundant AMR gene classes as being tetracyclines (35 ± 21 %), beta-lactams (28 ± 18 %), and MLS (27 % ± 17 %) ( Fig. 2 and S2B). Out of 238 AMR genes detected, the following 10 AMR genes were the most abundant at the shelter facility manhole (listed in decreasing order): tem, tetQ, tetW, msrE, tetO, ermF, tetX, mphD, aadA, and ermG.
Several important studies have emerged recently on how a large fraction of hospitalized COVID-19 patients during 2020 received antibiotics despite only a small fraction of admissions being associated with bacterial infections (Nori et al., 2021;Rawson et al., 2020;Lardaro et al., 2021). To test whether we could reach similar conclusions by the analysis of sewage and thereby develop an early warning system for AMR, we investigated AMR gene abundance levels in the WWTP during the COVID-19 surge in late 2020. Interestingly, we observed a spike in abundance (based on sequencing reads) for some AMR genes at the WWTP in mid-December 2020, followed by a short-term decline in late December 2020/early January 2021, and then another sudden but smaller spike in early to mid-January 2021 ( Fig. 3B and C). Interestingly, the corresponding wastewater SARS-CoV-2 viral count data and confirmed COVID-19 case counts also exhibited this characteristic 'double peak' during the winter surge in Las Vegas (Fig. S3).
Based on the timing of the surge in SARS-CoV-2 cases and how novel the disease was at that time, we asked whether the increase in abundance of AMR genes in wastewater was due to the administration of antibiotics associated with the treatment of upper-respiratory infections, which was an initial treatment option for SARS-CoV-2 infection (Excellence NIfHaC, 2020; Langford et al., 2021). Antibiotics associated with treating SARS-CoV-2 infections include beta-lactams (amoxicillin and ampicillin), MLS (azithromycin and clarithromycin), and fluoroquinolones (ciprofloxacin, moxifloxacin, levofloxacin) (Excellence NIfHaC, 2020; Langford et al., 2021;Rothe et al., 2021;Seaton et al., 2021). More AMR genes relating to beta-lactam (Fig. 3B) and fluoroquinolone (Fig. 3C) resistance were detected in the WWTP during the December 2020 spike in infections, while no significant difference was observed for MLS resistance genes (Table 1). Aminoglycosides and multi-drug resistance AMR gene classes also differed significantly over time, with aminoglycosides exhibiting significant changes from November to December and December to January (p-value < 0.05) and multi-drug resistance exhibiting significant changes from November to December (p-value < 0.05) (Table S3). Interestingly, the beta-lactam (Fig. 3E) and fluoroquinolone (Fig. 3F) AMR gene classes were more abundant in the shelter manhole in November 2020-at the onset of the COVID-19 outbreak within the facility (Fig. 3D)-compared to December 2020 when the larger surges in wastewater concentrations and confirmed cases were observed in the broader community ( Fig. 3A and Fig. S3). No samples were collected at the shelter facility prior to late November 2020 or after mid-December 2020 so it was not possible to directly compare data from the WWTPs and shelter facility across similar timeframes.

AMR genes relating to last-resort antibiotics
Carbapenems, vancomycin, and colistin are 'last-resort' antibiotics, which are used when first-line antibiotics are ineffective. These antibiotics are classified as 'last-resort' because resistance development and dissemination would make all current treatment options ineffective. Carbapenem resistance genes include ges-type, kpc-type, vim-type, ndm-type, tem-type, ctx-type and oxa-type carbapenemases (Du et al., 2018;Diene and Rolain, 2014;Poirel et al., 2019). Vancomycin resistance genes include vanA, vanB, vanC, vanD and several other genes that give rise to reduced susceptibilities (Arias and Murray, 2012;Ahmed and Baptiste, 2018). Colistin (polymyxin) resistance genes include the plasmid-mediated mobile colistin resistance gene (mcr), which is monitored closely because colistin is a 'last-resort' antibiotic used to treat carbapenem-resistant Enterobacterales (CRE) infections (Wang et al., 2018;Hussein et al., 2021). Many of these AMR genes were detected during this study, and more were identified in samples collected during mid-December from the WWTP (Fig. 4A) and late November from the shelter manhole (Fig. 4B). These differences were also statistically significant, suggesting that the abundance of 'last-resort' AMR gene classes was correlated with the surge in COVID-19 infections at the community level and the onset of the COVID-19 outbreak at the shelter. In addition, this was primarily reflected as an increase in the carbapenem resistance AMR gene class (Table 2).

A comparison of sewage AMR data with historical antibiograms from the region
To correlate our sewage data with the clinical documentation of antimicrobial susceptibilities for bacterial isolates, we analyzed community outpatient antibiograms for Las Vegas from 2018 to 2019. Interestingly, the most reported antibiotic class in the 2018-2019 antibiograms was betalactams, which accounted for 51 ± 3 % of reported antibiotics tested ( Table 3). The majority of Gram-negative isolates were resistant to common beta-lactams (amoxicillin and ampicillin), while a few were resistant to cephalosporins (cefazolin, cefoxitin, cefuroxime) (Tables S4-S7). Only E. faecium showed resistance to fluroquinolones, while other isolates showed mild susceptibilities. The antibiogram data indicate that basal levels of AMR genes relating to these classes should consistently be detected at local WWTPs, and intermittent spikes might reflect short-term changes in community conditions (e.g., increased prescription use during the winter 2020/2021 COVID-19 surge). While these antibiograms aid healthcare providers with decisions on empiric antimicrobial therapy, the significant  labor required to generate such datasets limits reporting to a small number of commonly isolated bacterial pathogens. In addition, the antibiograms report percentages of isolates susceptible/resistant to a particular antibiotic, and for the resistant isolates, the corresponding minimum inhibitory concentration is not determined. Similar clinical reports were not aggregated by the Southern Nevada Health District in 2020-2021 due to COVIDrelated logistical challenges. Therefore, it is not possible to directly compare the wastewater AMR profiles observed in this study with empirical culture-based clinical data across the same time periods.

Discussion
Studies have shown that 80-90 % of sewage-derived sequences originate from non-human origin and are unique to sewer infrastructure, but the remaining 10-20 % can be attributed to human fecal origin (Shanks et al., 2013;Vandewalle et al., 2012). Thus, data from sewage sequencing can theoretically be used to develop a community-level human gut microbiome profile. This approach has previously been used to assess whether a community is obese or lean, specifically by assessing relative abundance of Bacteroidaceae, Prevotellaceae, and Lachnospiraceae/ Ruminococcaceae (Newton et al., 2015). Bacteroidaceae and Prevotellaceae belong to the Bacteroidales order, while Lachnospiraceae/Ruminococcaceae belong to the Eubacteriales order. In the current study, Bacteroidales and Eubacteriales were the most abundant orders at the shelter (Fig. 1B) suggesting that human fecal microbiota may have dominated those facility-specific samples. Bacteroidales and Camplyobacterales were the most abundant at the WWTP (Fig. 1B) suggesting more diverse sources, which might be expected in the more complex sewage at a WWTP serving the broader community and traveling longer distances through the sewer collection system. Because of differing objectives, target-enrichment metagenomic sequencing is less discriminatory than 16S rRNA gene sequencing, but nonetheless, it was still possible to characterize microbial community structure in this study.
In contrast with a different study that identified macrolide resistant genes as the most abundant AMR class in North American sewage, our study observed macrolide resistance genes to be the second most abundant, with tetracycline resistance genes being the most abundant ( Fig. 2 and Fig. S2) (Hendriksen et al., 2019). A reason for this discrepancy could be due to the COVID-19 pandemic, or due to technical differences in the detection of AMR genes. For example, our extraction and sequencing approach simultaneously considers DNA (i.e., genes) and cDNA (expressed genes), while typical shotgun metagenomic studies focus only on DNA. In addition, the other study (Beran et al., n.d.) utilized a different method of detecting AMR genes where reads were mapped using MGMapper against the ResFinder database. As a result, the threshold for detecting AMR markers may contribute to how the data are ranked (Arango-Argoty et al., 2018;Ko et al., 2022).
Homelessness is a major societal concern, in part due to the potential for deterioration of urban and environmental quality coupled with poor public health conditions leading to high rates of disease transmission (Ly et al., 2019;Gerrity et al., 2022). Few studies have looked at AMR gene prevalence in this vulnerable population, but one study using nasal/pharyngeal samples reported lower abundance of AMR genes within the sheltered homeless population compared to controls (Ly et al., 2019). Sampling of homeless populations can be difficult due to the transient lifestyle and lack of participation, which constrains studies on the public health impact of homelessness. WBE conducted at shelter facilities may overcome these limitations by providing a more comprehensive and less invasive assessment of this disadvantaged, high-risk population. WBE could potentially target prevalence and transmission of viral infections (e.g., SARS-CoV-2, Fig. 4. Heatmap of AMR genes associated with carbapenem, colistin, and vancomycin resistance. AMR data from the (A) WWTP and (B) shelter were plotted over time to highlight the impact of the COVID-19 surge on 'last-resort' AMR gene classes. Genes were hierarchically clustered by abundances (cpm) using Pearson correlation coefficients. Colors represent natural log (ln) transformed abundances of different AMR genes.  hepatitis, enteric viruses) and other public health markers such as AMR.
Our results showed distinct peaks in AMR gene abundance at a WWTP and homeless/rehabilitation shelter, but the temporal patterns and profiles were different, perhaps due to differences in scale (i.e., community-level vs. facility-level) or demographics (i.e., broader community vs. sheltered homeless individuals). With respect to temporal differences, AMR gene classes associated with COVID-19 treatment exhibited greater abundance in the shelter sewage ( Fig. 3E and F) two weeks prior to those same trends being observed at the WWTP ( Fig. 3B and C). This suggests that community-level results may be useful for understanding public health conditions across the broader community, but more localized sampling may be warranted for targeted subpopulations. The SNHD antibiogram data (Tables S4-S7 and Table 3) did not provide actionable information because it only reported organisms with available testing data for ≥30 isolates and did not include organisms under surveillance, such as methicillin-resistant Staphylococcus aureus, Clostridioides difficile, and multi-drug resistance Neisseria gonorrhea. Furthermore, some of the antibiotics listed had no reported susceptibility percentages for many of the reported organisms. Acquiring local prescription data is more difficult for research purposes but is available to healthcare professionals, as well as in-depth antibiograms. WBE targeting AMR should not be considered a replacement for antibiograms but can complement more traditional culturebased AMR surveillance by providing more robust community-level AMR profiles through 'pooled' sampling, which might be more effective in alerting healthcare professionals of emerging threats. Recently, M. morganii isolated from a WWTP in Georgia, USA, was found to contain the plasmid-borne mcr-9.1 gene, which confers resistance to the lastresort antibiotic colistin (Hassan et al., 2021). This finding highlights the importance of sampling WWTPs for AMR genes and shows how this emerging AMR threat may be spreading in the U.S. This also demonstrates how wastewater surveillance may identify public health threats prior to their clinical confirmation, as was demonstrated previously for the detection of the Alpha variant of SARS-CoV-2 in Southern Nevada (Bolze et al., 2022).
This study included only a limited number of samples for generating the observed AMR profiles, so additional sampling is warranted to compare results in other locations and during other time periods of interest. Because of the widespread implementation of WBE for SARS-CoV-2, our experimental approach could be replicated for the same time period in other geographic locations using archived nucleic acid extracts. However, samples collected during subsequent COVID-19 surges (e.g., Delta and Omicron) may not complement our data due to the medical industry's evolving understanding of COVID-19 and approaches to treatment (e.g., less indiscriminate use of antibiotics). Nonetheless, WBE might be able to identify other temporal factors beyond COVID-19 affecting the efficacy of antibiotic stewardship programs. Thus, this study provides supporting evidence that WBE can provide valuable information about understanding and combating the proliferation of AMR.

Conclusion
This study demonstrates that a commercially available AMR target enrichment sequencing panel designed for clinical samples can be used to monitor AMR profiles in wastewater at the community and facility level. Increased abundance of AMR genes associated with beta-lactam and fluoroquinolone resistance was observed in wastewater samples collected during the winter 2020/2021 COVID-19 surge in Las Vegas, Nevada. Furthermore, this study demonstrates how wastewater-based epidemiology can be implemented to monitor genes associated with resistance to 'last-resort' antibiotics. In summary, this is a proof-of-principle study using a new NGS approach on non-clinical samples and provides a framework for future implementation of WBE to monitor emerging AMR threats.

Data availability
Data will be made available on request.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.