Wastewater surveillance reveals patterns of antibiotic resistance across the United States

Domain: biology
Published: 2026-04-01
Editorial outcome: major_revision

Abstract
Antibiotic resistance is a growing public health threat, with over 2.8 million antibiotic-resistant infections and 35,000 attributable deaths annually in the U.S. Here, we sought to use wastewater monitoring to assess community-level burden of antibiotic resistance. This study quantifies concentrations of antibiotic resistance genes (ARGs) by digital droplet PCR in wastewater solids obtained from 163 wastewater treatment plants in 40 states in the United States. We measure 11 ARGs that confer resistance to beta-lactams (bla , bla , bla , bla , bla , bla , bla , mecA), colistin (mcr-1), tetracycline (tetW), and vancomycin (vanA). The South has higher overall ARG concentrations compared to the Midwest. We pair these data with national data sets including antibiotic use, social vulnerability, size of animal agriculture operations, location of healthcare facilities, and presence of airports to investigate potential drivers of resistance. We also generate predictive maps of ARG concentrations for counties with data within the range of our training set in the United States. We show social vulnerability indicators (housing burden and access to health insurance) and indicators of international travel are associated with increased ARG concentrations in wastewater, while antibiotic usage is only weakly positively associated. Our results provide a national baseline of ARG concentrations and highlight the complexity of factors driving spread of antibiotic resistance.

Introduction
Antimicrobial resistance (AMR) is a critical and growing global public health threat, with resistant bacterial infections, such as carbapenem resistant Enterobacterales and vancomycin resistant Enterococcus, becoming more prevalent.1 The 2019 Antibiotic Resistance Threats Report by the U.S. Centers for Disease Control and Prevention (CDC) highlights the severity of the crisis, reporting over 2.8 million antibiotic resistant infections annually, resulting in more than 35,000 deaths.2 This is an underestimate, as it is based on people who seek medical attention.

Traditional clinical surveillance for antibiotic resistant infections involves testing bacteria isolated from clinical specimens to assess antibiotic susceptibility. This method is resource-intensive and only reflects antibiotic resistance genes (ARGs) in individuals seeking medical treatment, limiting our understanding of the true burden and diversity of AMR in the broader community.

Various anthropogenic factors contribute to the levels of AMR circulating in communities, but the relative influence of these factors remains unclear. Improper use of antibiotics in both healthcare and agriculture has long been recognized as a major driver of antibiotic resistance, with overprescription and misuse accelerating the development of resistant bacteria.3 However, recent studies have increasingly emphasized the role of socioeconomic factors in shaping resistance patterns.4,5 Disparities in access to healthcare, sanitation, and education can influence not only antibiotic use but also the effectiveness of infection control measures. In a global analysis, Collignon et al. reported limited association between clinical isolates, antimicrobial resistance levels, and antibiotic consumption but found stronger association with socioeconomic factors such as GDP per capita, education, infrastructure, and public healthcare spending.5 Similarly, a global analysis of sewage metagenomes by Gupta et al. and Lee et al. identified significant correlations between resistance levels and socioeconomic factors.6,7 Additionally, global travel has also been found to play a key role in spreading resistant pathogens, facilitating the dissemination of resistance across borders.8,9 These interconnected drivers highlight the need for a comprehensive examination of multiple determinants of antibiotic resistance.

Recent advances in wastewater monitoring for SARS-CoV-2 and other viral pathogens have established data collection infrastructure that can be leveraged to monitor various markers of infectious disease and human health, including ARGs, in a timely and comprehensive manner.10 Wastewater has consistently been recognized as a reservoir for ARGs. While previous research has characterized ARG levels and their fate within wastewater treatment systems using PCR-based methods, studies undertaking national-scale community antibiotic resistance surveillance via wastewater monitoring are still emerging. Prior work assessing community-level ARGs in wastewater has largely involved either global spread or differences in the wastewater resistome globally,4,6,11,12 or are limited to single cities or states.13–15 Additionally, many previous wastewater ARG studies conducted for the purpose of monitoring community antibiotic resistance have used shotgun metagenomics. While metagenomic approaches are useful in understanding the overall ARG diversity, they often lack the sensitivity and quantitative precision needed for accurate abundance estimates. In contrast, digital PCR (dPCR) has the advantage of being highly sensitive and quantitative, allowing for normalized abundance estimates and comparison across samples in time and space.16–18

In this study, we use droplet digital PCR (ddPCR) to quantitatively measure the abundance of 11 clinically relevant ARGs that confer resistance to beta-lactams (bla , bla , bla , bla , bla , bla , bla , mecA), colistin (mcr-1), tetracycline (tetW), and vancomycin (vanA) from 163 wastewater sites to generate a cross-sectional data set of the U.S. By linking mostly publicly available secondary data, we investigate potential determinants of wastewater ARG concentrations, including antibiotic use, social vulnerability, the size of animal agriculture operations, and the locations of healthcare facilities and airports. Finally, using a subset of our data as a training set, we construct a random forest model to predict antibiotic resistance gene prevalence in counties across the U.S. that can be interpolated with our training set.

Results
Clinically relevant ARGs are at high concentrations in wastewater across the U.S.

We measured the concentration of 11 antimicrobial resistance genes (ARGs) across 163 sites in the U.S, spanning 40 states for a cross-sectional data set (Figure 1A). An average of 2.7 samples were collected per wastewater treatment plant (WWTP) over a week in May of 2024 and concentrations were averaged to obtain one representative concentration for each target per treatment plant (Table S1). We selected the targets based on the Centers for Disease Control Antibiotic Resistance Threats Report that identifies pathogen-resistance combinations as urgent, serious, and concerning based on their threat to human health.1 Targets included genes that confer resistance to beta-lactams, colistin, tetracycline, and vancomycin. The ARGs bla , bla , bla , bla , bla , and tetW consistently exceeded 10-7 copies per copy of 16S rRNA gene and were detected in all samples (Figure 1B). Among these, tetW exhibited the highest concentration across all samples, with minimal variation in its distribution. In contrast, bla , bla , mcr-1, mecA, and vanA displayed more variability in their distribution. Notably, bla and bla concentrations exhibited a broad distribution, ranging from undetected to 3x10-5 and 8x10-4 copies per 16S rRNA gene respectively. All targets in this group showed either non-detects (i.e. all wastewater samples collected during the week of sampling were non-detects for the WWTP) or partial detects (i.e., at least one sample collected for the WWTP during week of sampling was a non-detect). Mcr-1 had the highest number of non-detects (38.3% of sites) and partial detects (32% of sites), followed by vanA (24.5% of sites non-detects and 22.7% of sites partial detects) and bla (21.5% of sites non-detects and 16.6% of sites partial detects). Overall, most 16S rRNA normalized ARGs in wastewater were positively correlated with each other (Spearman’s correlation coefficients ranged from 0.17 to 0.88, all p < 0.05), with beta-lactamase genes showing some of the strongest pairwise correlations, except for bla (Figure S1A).

To summarize the resistance burden, we calculated z-scores based on the distributions of each gene normalized by 16S rRNA gene concentrations, generating burden scores for each resistance class and total resistance. Although more than the presence of ARGs is needed to estimate burden, we will use this term here to mean antibiotic resistance gene prevalence in the community’s wastewater. Burden z-scores were highly correlated (Figure S1B). Comparison of these burden scores across census regions revealed significant regional difference in total resistance burden (Kruskal-Wallis test, p < 0.05; Figure 1C). Specifically, the South had a significantly higher total resistance burden compared to the Midwest. Additionally, the West and the South demonstrated significantly higher colistin resistance burdens than the Midwest and the Northeast (Figure 1D; Post-hoc Conover-Iman test, p < 0.05, Benjamini-Hochberg correction applied). For individual beta-lactamase genes, significant regional differences were observed without a clear trend (Figure 1D). For instance, the Midwest had significantly higher bla levels but lower bla concentrations (Figure S2). Regional differences among all ARGs can be seen in Figure S2. The two genes with the highest national variability, bla and bla , displayed distinct patterns: bla showed significant regional differences, with the South and the West having higher concentrations than the Midwest, while bla exhibited similar variability across all regions.

Determinants of wastewater antibiotic resistance burden

Using publicly available datasets, we investigated the influence of specific anthropogenic or socioeconomic factors on the antibiotic resistance burden observed in wastewater. Values for secondary variables for sites were determined at either the sewershed scale or for the county predominantly intersecting a sewershed (see Supporting Information). From our correlational analysis between resistance gene prevalence and potential determinants, several notable associations were identified. Colistin resistance measured in wastewater was associated with overcrowding, speaking limited English, and certain race and ethnicities (p < 0.05 after Benjamini-Hochberg adjustment; Figure 2A). Overall beta-lactamase resistance was associated weakly with many variables but specific beta-lactamase genes (bla , bla , bla , bla ) showed stronger correlations with socioeconomic factors such as being uninsured, being under housing burden, overcrowding, and lack of a high school diploma, as well as speaking limited English and race and ethnicity variables (p < 0.05 after Benjamini-Hochberg adjustment; Figure 2B). Among these four genes, bla exhibited the greatest number of significant correlations with socioeconomic factors. Antibiotic use from Epic Cosmos and the number of animals had little to no correlation with the resistance burden. These correlations were adjusted with Benjamini-Hochberg method to control the false discovery rate (Figure 2); unadjusted correlation coefficients are provided in Figure S3.

Total antibiotic resistance burden score differed significantly based on the presence of hospitals and whether the number of nursing homes in the sewershed exceeded the median (4) (Figure 2C). While colistin resistance gene concentration normalized by 16S rRNA did not differ based on the presence of hospitals or nursing homes, higher mcr-1 concentrations were associated with the presence of airports, higher population density, and urbanicity (Figure 2D). When examining individual beta-lactamase genes, some genes, such as bla showed no significant difference with presence of points of interest. In contrast, other genes like bla , bla , bla , and bla showed significant differences based on presence of hospitals and the number of nursing homes (Figure 2E, F; Figure S4E, H, I). Notably, higher bla concentration was the only beta-lactamase gene that was associated with the presence of airports, higher population density, and urbanicity (Figure 2F).

Predictors of antimicrobial resistance concentrations

Using publicly available datasets of non-clinical determinants of antibiotic resistance and antibiotic prescription rates obtained from EPIC Cosmos, we developed random forest models to predict wastewater ARG concentrations across U.S. counties with characteristics similar to our sampled wastewater sewersheds. The random forest models presented here had one tuning parameter: the number of variables randomly sampled at each split (node), which ranged from 2 to 28. For all ARGs, the model that minimized RMSE considered 2 variables at each split, so all subsequent results are presented with this configuration. Overall, the predictive performance for individual genes, indicated by R2, ranged from 0.14 to 0.46 (Figure S5). The mcr-1 model, showing colistin resistance, achieved one of the highest predictive performances with an R² of 0.46, indicating it could explain 46% of the variance in colistin resistance across counties based solely on the secondary dataset. The models for bla and bla , both conferring carbapenem resistance and classified as urgent threats by the CDC, showed similar predictive performance with R² values of 0.30 and 0.37 respectively. The vanA model, showing resistance to vancomycin, had an R2 of 0.34 (Figure 3). Variables identified to be the most important across all individual ARG prediction models were proportion of the population being under housing burden, being uninsured, certain race and ethnicities, speaking limited English, and urbanicity of the sewershed (Figure 4).

Methods
The authors confirm that this research complies with all relevant ethical regulations. Wastewater sampling and analysis did not involve human participants or the collection of identifiable private information; therefore institutional review board oversight was not required for this component of the study. Access to antibiotic prescription data used in this study was obtained through the Epic Cosmos platform. In accordance with the Epic Cosmos data use and publication policies, the study protocol followed the established approval process, which included formal notification of the intent to publish and the submission of the manuscript for review and feedback prior to submission.

ddPCR assay optimization

We chose probe-based qPCR assays in literature that have been utilized for detection of the target ARGs in wastewater to test on the ddPCR platform (Supplementary Data 1). The assays were initially tested on gblocks (IDT; Supplementary Data 1) of the target region of each ARG as either a singleplex (FAM) or a duplex (HEX) assay with an annealing temperature gradient ranging from 55°C to 65°C in order to identify an annealing temperature that yielded sufficient separation between negative and positive droplets (Figure S9). Subsequently, the assays were tested using templates at four different dilutions (neat, diluted in molecular grade water at 1:10, 1:100, 1:1000) of wastewater samples collected from seven wastewater treatment plants in the U.S. with distinct geographical characteristics to identify the dilution factor suitable for each target ARG. Consequently, we categorized the target ARGs into two groups based on their expected concentrations in wastewater (high versus low) in order to multiplex the assays into two 6-plex assay mixes that utilized an annealing temperature of 58°C.

Sample collection

Wastewater treatment plant (WWTP) staff provided either “grab” samples from the primary clarifier or 24-hour composite influent samples from the headworks during a one-week period in May 2024 (Supplementary Data 2 shows which WWTP collected each type and the date range of the collected samples). The “grab” samples in this case are solids collected in the primary clarifier over 1-8 hours, representing composite community wastewater solids. The samples were then stored at 4oC, shipped to the laboratory overnight, and processed within 48 hours of receipt at the laboratory. A total of 443 samples were collected and analyzed for this study (Figure 1A, Table S1). The total population served by the treatment plants was 24,455,284, representing 7% of the US population.

Pre-analytical processing and nucleic-acid extraction

Our study focused on the analysis of wastewater solids for ARG quantification. This decision was informed by previous comparative studies, including our own work, which have demonstrated that ARGs are reliably detectable and highly correlated across different wastewater matrices, including raw influent, centrifuged influent pellets, and wastewater solids. Specifically, our prior research showed comparable detection sensitivity and strong correlation in ARG concentration between these sample types, with wastewater solids exhibiting slightly higher ARG concentrations on a per-mass basis.44 This indicates that wastewater solids serve as a representative and effective matrix for assessing community-level ARG prevalence and that centrifuged influent pellets were a suitable substitute when wastewater solids were not available.

Methods for pre-analytical processing and nucleic-acid extraction and purification have been described in detail by Boehm et al., so are briefly summarized here.45 Wastewater solids were obtained from the sample by use of Imhof cones for influent samples and centrifugation for primary settled solids and then suspended in a buffer (DNA/RNA shield, Zymo Research, R1100-250, California, USA) containing bovine coronavirus (as a positive extraction control; Calf-Guard Cattle Vaccine, Ohio, USA) at a concentration known to reduce the potential for inhibition in the analytical steps.46 The samples were homogenized after addition of grinding balls (5/32-in., OPS Diagnostics, New Jersey, USA) and then nucleic-acids extracted and purified from 6-10 aliquots of the supernatant to obtain 6-10 300 µl replicate extracts from each sample (see Supplementary Data 2 for whether 6 or 10 replicates used for the plants). Nucleic-acid extraction and inhibitor removal was conducted using Chemagic Viral DNA/RNA 300 kit H96 (Revvity, CMG-1033-S, Massachusetts, USA) for the Chemagic 360 and Zymo OneStep-96 PCR Inhibitor Removal kit (Zymo Research, D6030, California, USA). The methods for nucleic acid extraction are described step-by-step in protocols.io.47 The first two aliquots of every sample were chosen to be analyzed for ARGs. Extraction-negative controls (water) were extracted using the same protocol as the homogenized samples. Nucleic acids were used immediately for quantification of ARGs. A portion of the solids were dried to determine their dry weight.45

Droplet Generation and Quantification

The probe used for each assay was labeled with one of six fluorophores (FAM, HEX, Cy5, Cy5.5, ROX, and/or ATTO590) and run in a 6-plex assay using the probe-mixing approach (Supplementary Data 1). All targets were multiplexed and run together based on being high concentration (samples diluted 1:100 with ultrapure H O) or low concentration (undiluted samples) in wastewater. Two wells were run per reaction using the first two nucleic-acid extractions. Digital droplet PCR was performed on 20 µL samples from a 22 µL reaction volume, which consisted of a 5.5 µL template. Final primer concentrations were 900 nM and probe concentrations were 250 nM for all ARG targets. All ARG master mixes were prepared using ddPCR Supermix for Probes (BioRad #1863010, California, USA). Droplet generation was performed using an Automated Droplet Generator (BioRad #1864101, California, USA). Quantification was performed using a QX600 Droplet Digital PCR system (BioRad #17007769, California, USA). All assays were performed at an annealing temperature of 58°C. Thermocycling conditions were as follows: 5 minutes at 95°C, followed by 40 cycles of 30 seconds at 95°C and 1 minute at 58°C, followed by 10 minutes at 98°C, and an infinite hold at 4°C. Each plate run included no template controls, extraction negative controls, extraction positive controls, and positive controls, which all performed as expected. Extraction negative controls used ultra-pure water, extraction positive controls used BCoV spiked into each wastewater sample, and positive controls were synthetic IDT gblocks which included the target amplicon. No template controls and extraction negative controls were negative for all targets in all runs. Extraction positive controls and template positive controls were positive in all runs. For each well in ddPCR, a minimum of 10,000 droplets was required and at least three positive droplets were required to be quantifiable. If a target was undetected from a diluted sample, the measurement was redone or excluded (one mecA, six bla OXA-48 , and one mcr-1 measurement excluded). In cases of duplicate measurements at different dilutions, the value with the smaller ddPCR confidence interval was selected. Two replicate wells were merged per measurement. Thresholding for each target was performed automatically using the QX Manager Software Version 2.2 (BioRad, #168032, California, USA).

After quantification of ARGs, the leftover nucleic acid replicates were combined and stored at -80°C for 17 days to 21 days prior to being used as templates for measuring high copy number bacterial target, 16S rRNA. 16S rRNA genes were quantified using ddPCR EvaGreen Supermix (BioRad, #1864034, California, USA) and QX200 Droplet Digital PCR system (BioRad, #1864001, California, USA). EvaGreen is an intercalating dye that binds to double stranded DNA. Final primer concentrations were 100 nM for 22 µL reaction volume, 5 µL of which was the template. 16S rRNA quantification was run in singleplex reactions using template diluted 1:50,000 with UltraPure H O. One well was run per reaction using the combined template from the quantification of ARG genes. Thermocycling conditions were as follows: 5 minutes at 95°C, followed by 40 cycles of 30 seconds at 95°C and 1 minute at 55°C, followed by 5 minutes at 4°C, 5 minutes at 90°C and infinite hold at 4°C. Each plate run included no template controls (NTCs). The NTCs resulted in a small number of positive droplets (average 16 positive droplets) and therefore any sample concentration on the same plate with less than an order of magnitude difference was excluded from analysis (three 16S rRNA measurements excluded). All wells had more than 10,000 partitions. Thresholding for each target was performed automatically using QX Manager Software Version 2.2.

Additional details associated with the Minimal Information MIQE Experiments (dMIQE) for digital droplet PCR reporting are described here.48 The average (standard deviation) number of droplets for the quantification of ARGs was 34225 (7143; 2 merged wells) and for bacterial count was 17997 (1443; 1 well). The volume of partitions, as reported by BioRad, was 0.000795 µL. The mean and standard deviation of copies per partitions for each target is shown in Table S4. BCoV served as an extraction and inhibition control to assess nucleic acid recovery efficiency and to identify samples with significant PCR inhibition. Median bovine coronavirus (BCoV) recoveries (BCoV measured in sample / BCoV measured in DNA/RNA shield without sample) were, on average 1.06 (25th percentile = 0.71, 75th percentile = 1.3) across all samples. Recovery higher than 100% is likely due to an underestimation of the amount of BCoV spiked into the samples. Based on the recovery, we concluded that there was no gross inhibition of the samples and no samples were removed due to low BCoV recovery. After QA/QC, we had 443 different samples and 5309 measurements that represented 163 WWTPs in 40 states and 134 cities in the U.S.

Demographic, socioeconomic, and environmental data

Antibiotic resistance may be influenced by various demographic, socioeconomic, and environmental determinants—in addition to antibiotic usage. Demographic and socioeconomic factors affect a community’s capability to respond and adapt to external hazards and stressors, such as public health emergencies. A community’s status as urban or rural or its population density may affect its access to and quality of healthcare services. Proximity to major points of entry (e.g., airports) and health care facilities (e.g., hospitals, nursing homes) may facilitate the importation or transmission of public health threats in communities. Antibiotics are also commonly used for agricultural purposes, and ARGs may be transmitted zoonotically within communities. Therefore, we characterized sewersheds based on several demographic, socioeconomic (Table 1), environmental, and antibiotic usage determinants (Table 2) to assess their relationship with wastewater concentrations of ARGs.

Most WWTPs (n = 94 of 163) provided a sewershed boundary. For WWTPs that did not provide a sewershed boundary (n = 69 of 163), we approximated a boundary based on the zip code(s) serviced by the WWTP and 2023 USA ZIP Code Boundaries published by Esri (Source: TomTom, US Postal Service, Esri).49 Because sewershed boundaries do not align with traditional boundaries (e.g., counties or census tracts), we approximated sewershed-level estimates for each determinant in Table 1 and 2 as described in the Supporting Information. The demographic and socioeconomic determinants we selected from the American Community Survey (ACS) are based on the 2022 Social Vulnerability Index (SVI) from the US CDC as social vulnerability describes the degree to which demographic and socioeconomic factors affect a community’s resilience.50,51 We derived urbanicity and population density from the US Census Bureau and ACS, respectively.51,52 Sewersheds with an urban area proportion ≥ 50% were classified as urban and sewersheds with an urban area proportion < 50% were classified as non-urban for data analysis. We obtained locations of major airports from Esri (source: Federal Aviation Administration’s National Airspace System Resource Aeronautical Data Product) and locations of hospitals and nursing homes from the US Department of Homeland Security’s Homeland Infrastructure Foundation-Level Data database.53,54 We obtained animal inventory numbers from the US Department of Agriculture National Agricultural Statistics Service.55 Antibiotic prescription was determined using Epic Cosmos (Epic Cosmos, Epic Systems Corporation, Wisconsin), a dataset created in collaboration with a community of Epic health systems representing more than 284 million patient records from over 1,500 hospitals and 36,000 clinics from all 50 states, D.C., Lebanon, and Saudi Arabia, though only U.S. data was used in this study. .56,57 The Epic Cosmos data includes all healthcare systems which use the Epic electronic healthcare record system. This has been shown to be nationally representative from a demographics perspective and nationally representative of tested sub-populations.57,58 To calculate antibiotic usage, all encounters—which represent any visit, outpatient or inpatient, between a patient and healthcare provider—were selected for every county available in Epic Cosmos between June 1, 2023 and June 1, 2024. From the selected encounters, the number of prescriptions for each antibiotic class was calculated for each county, and the total number of encounters in the county were used to normalize the data. A total of 1,590,238,036 patient encounters were analyzed.

Data analysis

Analysis was conducted using R (Version 4.3.0), ArcGIS Pro (version 3.1.1), and the following packages: caret (version 7.0-1), dplyr (version 1.1.2), EnvStats (version 2.8.1), corrplot (version 0.92), conover.test (version 1.1.6), stats (version 4.3.0), ggplot2 (version 3.5.1), and pheatmap (1.0.12). Outlier detection was performed using Rosner’s test on 16S rRNA to exclude significant outliers. One 16S rRNA measurement was excluded. Each time point was treated as a replicate for each wastewater treatment plant and therefore averaged to obtain a final representative concentration of each target. This approach was specifically designed to mitigate the influence of short-term fluctuations in wastewater characteristics, thereby providing a pragmatic method for a more robust and representative average concentration for each site during the cross-sectional sampling period. Samples were categorized as non-detect (ND) if none of the replicates detected the target, detect (D) if all replicates detected it, and partial detect (PD) if at least one of the replicates did not detect the target. For non-detect values, half of the lower measurement limit was calculated and assigned as a substitute value.59 Summary of the analysis when 0 was substituted for non-detect values is shown in Figure S10. All ARG concentrations were normalized by 16S rRNA concentration to account for variations in bacterial biomass in each sample.60,61 Z-scores were calculated for each target by comparing the concentration to the overall distribution across all samples, using the mean (μ) and standard deviation (σ) in the following equation: z = (x - μ) / σ. The AMR burden was calculated by summing up the z-scores of all ARG targets. The beta-lactamase burden was computed by summing up the z-scores of bla , bla , bla , bla , bla , bla , and bla . Data points with z-scores greater than 3 or less than -3 were excluded from analysis involving secondary datasets (four bla , five bla , three bla , two bla , two bla , two bla , three bla , six mcr-1, one mecA, one tetW, and two vanA measurements excluded).

For statistical analysis, the Shapiro-Wilk test was used to assess normality, and Spearman’s correlation coefficient was employed to evaluate relationships between variables, as the data was not normally distributed. For determinants of resistance with insufficient distribution (number of airports, number of hospitals, number of nursing homes, population density, and urbanicity), the data was split into two groups to conduct a bivariate analysis using the Wilcoxon rank sum test to determine if there was a significant difference in wastewater ARG concentrations between the lower and upper values of potential determinants. Number of airports and hospitals were split into groups based on the presence of the facilities, and the number of nursing homes were split by the national median number of nursing homes (median: 4). The sewersheds were considered urban if there was 50% or more area considered to be urban. Population density was split by 617 per km2, which is the average number of people in cities in the U.S. The number of animals was analyzed as percentile ranking from the national dataset. Principal component analysis was conducted on determinants of resistance to visualize covariance among the secondary variables. The Kruskal-Wallis test was applied to assess differences across Census regions, followed by the Conover-Iman test for post-hoc analysis to identify significant pairwise differences. The Benjamini-Hochberg correction was applied when potential determinants were involved and when investigating regional differences to control for the false discovery rate.

Predictive modeling

Based on our secondary data’s characteristics, we chose random forest modeling to predict wastewater ARG concentrations in every county in the United States. Random forest modeling allows us to use data with varied distributions and nonlinear relationships, is resistant to outliers, provides interpretable feature importance, and is robust to covariance between independent variables.62,63 Given our limited sample size for the cross-sectional analysis, the random forest model minimizes the risk of overfitting often associated with more complex architectures (such as neural networks and gradient boost models) when trained on small datasets. We first aggregated our mean wastewater ARG measurements at the county level. If a county was represented more than once in the dataset (i.e. contained 2 treatment plants), one of the sets of measurements from that county was randomly selected to represent the entire county. We chose this approach instead of averaging the values to minimize structural bias, as averaging would systematically and incorrectly assign equal weight to all plants regardless of their population size, which we could not correct for since wastewater treatment plants typically do not report the number of population served broken down by county boundaries. A sensitivity analysis regarding this approach can be found in the SI. This left 108 unique observations. All variables from Table 1, excluding number of airports, were considered and obtained for each county in the U.S. Airport data was omitted from this analysis due to a highly skewed distribution in which almost all included sewersheds had no airports (93% of training data locations, 97.5% of nationwide locations). The variables for number of nursing homes, number of hospitals, and number of animals were turned into a percentile rank before analysis. To train the models, we used the caret package (version 7.0-1) in R.64 To avoid extrapolating extensively to rural or unsewered areas, we limited the predictive model to counties that have secondary data values within our minimum and maximum ranges of the training data set for all variables. For example, if our training dataset included 0 to 1000 nursing homes, only counties that had 0 to 1000 nursing homes were included. Ten-fold cross-validation was used to evaluate the models, and we varied the number of features at each split between 2 and 28. Models were selected based on minimizing root mean square error. Variable importance was assessed using the varImp() function from the caret package, which randomly permutes variables and assesses the impact on prediction accuracy (Figure S6). More details on our modeling approach can be found in the SI.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The contribution by Kim et al. has many strengths, but also some considerable weaknesses. The strength of the work lies in their extensive dataset of ARG concentrations from sites located across the US. The authors are clearly able to show that there is substantial variability in the concentrations of ARGs in different regions of the US. This is highly useful information that can help inform our collective understanding of ARG dissemination. The weaknesses of the work result primarily from the authors' attempt to not only consider the local sewershed-scale determinants of AMR, but to then extend their results nationally. As noted in my general comments #2-4 (below) there are many questions raised by their approaches to define the variables included in their assessments of the local importance of the different determinants. The US-scale simulations of ARG concentrations that are central to this work are fully dependent upon the quality of the data used to develop their random forest models and as such that data and its interpretation must be fully transparent.\nAt the present time, there appears to be too much extrapolation of what they have learned from their 163 sites (which are geographically skewed towards a few specific regions (primarily California)) to the much broader US. Relatedly, a further challenge to the potential validity of their US assessment is that their dataset is restricted to regions served by sewer networks; however, substantial portions of the US are unsewered and as such the relative importance of such unsewered sites (which often serve lower income regions) cannot be estimated using sewer-based data alone. These potential limitations (particularly given the relatively low R^2 values reported for their models) require further consideration. One potential approach might be for the authors to more geographically limit the modeling to specific regions characterized by large numbers of sampling sites.\nGeneral Comments:\n1. Much of the work described in this contribution is reliant upon the development of correlations between measured ARG levels and metadata obtained from outside sources. At the end of the day, even if the correlations are strong, they remain correlations and cannot be used to infer causation. Throughout the manuscript the authors need to include statements that reflect this fact. Further, when discussing these correlations the manuscript needs to better incorporate consideration of co-correlation in cases where the different metadata parameters are subject to significant covariance.\n2. Additional information needs to be provided to better explain many of the variables included in Table 1. How were the numbers of airports, hospitals, and nursing homes considered? Was there a local regional area considered or do they have to lie within the expected sewershed? It seems straightforward to consider the hospitals and nursing homes within the geographical area that overlays the sewershed; however, the number of airports would seem to be less easily defined since a sewershed may not overlap with an airport even though the people within the sewershed heavily use a nearby airport. Further, how is the size of the airport, hospital, or nursing home taken into account? One would anticipate that larger facilities might have a much more substantial impact on ARG levels detected in sewage.\n3. Based upon the data presented in Figures S12-S13 it appears that a majority of the sewersheds being considered contain no airports, no hospitals, no nursing homes, and no cattle, chickens, or hogs. Further, in the case of the three agricultural animals only a limited number of sites contain either chickens or hogs, but with substantial variation in total numbers (i.e., between 10^7 and 0) . Given the relative paucity of sewersheds with any of these variables it seems very unlikely that their evaluation of their role as ‘determinants of wastewater antibiotic resistance burden’ [Lines 158-186] can be accurately assessed and then accounted for in their subsequent US model development. Nonetheless, the authors attempt to do so as indicated by their text on lines 178-180: ‘…higher mcr-1 concentrations were associated with the presence of airports…’ and lines 269-271 'We also identified international travel and recent immigration...as significant drivers of increased ARG burden. NDM and colistin resistance were significantly associated with the number of airports...' [As an aside - how did the authors determine that international travel and immigration were important? Solely based on nearby airports? If so, that would seem to be quite an intellectual leap.] Based upon Figure S12 only ~5 sites had any airports at all. It would seem more appropriate to not consider these variables when looking at potential associations with ARG levels (or total burdens). Similarly, by extension, the US modeling should be re-evaluated without consideration of these variables.\n4. One of the potential key contributions of this manuscript is the modeled prediction of ARG prevalence across all counties in the U.S. However, it is not clear exactly how this model was constructed given the limited size of the initial dataset (163 wastewater sites with a majority present within California). As the authors illustrate throughout the manuscript, the factors that drive the dissemination of antibiotic resistance are highly complex, often quite inter-related, and site-specific. Accordingly, it seems to be quite the intellectual leap to think that their model is able to accurately predict ARG levels for regions of the country that are not included within their model training set. As illustrated in Figures S10 and S11 the distributions are heavily skewed towards urban areas that are only moderately populated – thus, the capability of their model to properly predict what is happening in either rural (potentially unsewered) or dense-urban cities would seem to be low.\n5. One of the key sources for metadata is the Epic Cosmos database. Unfortunately, the authors description of how this database was used is insufficient. Epic Cosmos is a semi-closed resource that not everyone (including this reviewer) has access to and thus it is unclear what types of data it provides. The SI makes use of the term ‘encounter’, but this term is not defined. Is it a drug prescription? Or something else?\nAdditional Specific Comments:\n1. Line 27. Is it necessary to make the statement starting with ‘This is an underestimate…’? Such statements may belong in the main body of the paper, but it is not clear it belongs in the abstract which should focus on the results produced as part of the work.\n2. Lines 66-67. Include commas after ‘clinical isolates’ and ‘antimicrobial resistance levels’.\n3. Lines 78-79. The the use of the word ‘most’ ‘…most research to date has focused on the fate of ARGs across the wastewater treatment processes’ is incorrect. There is a large body of research that has reported ARG levels in raw sewage over many years. Similarly, the statement ‘have mostly been conducted using shotgun metagenomics’ is also incorrect. There have been numerous reports wherein PCR-based methods have been used to quantify ARGs in raw sewage.\n4. Line 85. The authors suggest that metagenomics cannot detect ‘rare targets’ but provide no literature justification of that point. The typical assumption is that PCR provides lower levels of detection, but since it is a targeted approach it is unable to detect untargeted analytes. I suggest rewriting to highlight the better detection limit for PCR-based methods and not make any claims about ‘rare targets’ – particularly since such a use of PCR is not revisited later in the manuscript.\n5. Lines 101-123. It would be useful to compare the results presented in this manuscript to those obtained previously. How do the ARG levels compare? Should the reader expect that CMY, CTX-M, KPC, OXA-48, TEM, and tetW are ubiquitous? Should tetW be present at the highest level?\n6. Line 125. Not all readers may be familiar with ‘z-scores’ and thus it would be useful either here or in the Materials and Methods to mathematically define how they are calculated.\n7. Lines 176-186. How is urbanicity defined? The SI considers the distribution of urbanicity along a scale from 0 to 100, but Figure 2 defines urbanicity using a Yes/No scale. Is >50% urban a Yes or is something closer to 90% considered urban? The clear definition needs to be provided. A similar lack of clarity also affects the authors use of the term ‘urban’ – Figure 2 defines it using a Yes/No scale, but what that means in the context of Figure S11 is unclear. How dense does an area need to be to be considered urban?", "reviewer_2": "Thanks for asking me to review this manuscript. It was a pleasure to read.\nIn brief, Pickering et al., sought to understand factors that may contribute (associate with) to the distribution of AMR (in particular 11 targeted ARGs) across the US by measuring their abundance at 163 separate WWTP in the US. These sites included those within 40 separate states. They integrated this data with publicly available data on antimicrobial use and census data to look in a holistic fashion on factors that associate with AMR.\nThe manuscript is particularly well written and easy to follow. The use of the same colour for different regions made the manuscript easy to read. Its broad scope and relevance make it an ideal candidate for consideration of publication in Nature Communications. Perhaps the most contentious (only in some circles) is the finding that local antimicrobial prescribing (in this regulated, controlled environment) only have a modest impact. This is really important and a bit more discussion on this (what is there is great) would be fantastic given the magnitude of this global crisis.\nThere are several issues that should be addressed to make the manuscript more robust and easier to follow.\nMajor issues.\n1. Details on the sewersheds that are being monitored and when they are being assessed are key – and are lacking here.\na. Please include a table on the nature of each of the WWTP - and the size of the population served. As the authors rightfully highlight the comprehensive nature of their study (i.e. 40 states included!) - it would be nice to see what fraction of the US population is accounted for in these samples which would allow the reader to determine how much stock to put into their projections across each county.\nb. What are the influences of rainwater in the system? Some communities have separate storm and sewage wastewater systems - and although many do not. While 16S will somewhat mitigate this - it should be commented on.\ni. The authors use the term grab samples – but also describe this as potentially being an aggregate sample collected over 1-8. Please distinguish if this was the same or different – and include in the table relative to those that provided 24-hours composite samples. Are differences observed between them based on diurnal patterns of water usage?\nc. Samples appear to have been connected at a single point (or week in time) - are these the same through the entire study? If not - please indicate when for each site. What are there seasonal influences?\ni. While authors have argued some studies have not shown seasonal influence of AMR distribution – many, many have particularly as it pertains to respiratory viral seasons influence\nii. Limited serial monitoring with only 2.7 samples per WWTP all collected (I think) within one week. The authors do not describe or include information as to whether this was done at the same time and what time this might be. While some m\nd. Social variables in these WWTP require more details.\ni. Airports - what was defined as an airport? Are these only large hubs - by some specific definition or does it include even small non-commercial strips? While definitions report # of hospitals the author’s analysis is presented as presence/abscence.\nii. Hospitals/LTC - how are hospitals defined? Are they only tertiary or quaternary sites or specifically large LTC facilities? I find it hard to imagine that any of these large catchments for WWTP don’t include a small community hospital or nursing home? While these are presented as < or > 7 – what does this mean? 1 large nursing home of 350 individuals is likely more relevant than 8 small nursing homes of 40 individuals….\niii. More detail on all of these variables would help the reader.\n2. MecA is not - as stated - a beta-lactamase and should not be grouped with others. If the authors wish to continue to include it – they will need to modify several figures. Its inclusion, however, in this list as an AMR threat is a bit misleading. MecA is relevant in so far as its presence in S. aureus is a threat, however, the vast majority of its measured abundance will derive from CONS which would not be considered an AMR threat.\n3. Quality control is important – but how relevant is supplemented BCoV - an RNA target well established for SARS-CoV-2 - as a control in an AMR experiment? No results are presented on this – and I suggest to remove it entirely as it does not seem to be relevant at all. If you think so – please describe why it is important and provide data demonstrating which samples were excluded owing to the results generated.\n4. I am unfamiliar about EPIC Cosmos and its ability to provide comprehensive antimicrobial usage data. Is this human only - or does it include animals and agriculture? Does it include hospitals or just communities? Is it represent ALL antimicrobial prescribing or just a small sampling? - please provide some background even if it is in the supplement for the reader to understand. Currently supporting citation is the company.\na. From this data the authors have assessed beta-lactam usage. However, lumping all beta-lactam antibiotics together for assessing selection influence for AMR is likely inappropriate. Beta-lactams can be very narrow- or very broad- spectrum. What the authors have done is like lumping all things with wheels together …. Some with wheels might be unicycles others semi-trailer trucks and others yet are 747s…. It is not the wheel (i.e. the beta-lactam ring) that differentiates them. The authors should at the very least look at narrow (e.g. penicillin, amoxicillin, ampicillin) vs broad (amoxi-clav, amp-clav, pip-tzp, ert, mem, imi) etc separately.\n5. The authors have identified the highest burden of AMR is in the NE and the S – and this is clear in Figure 1. However, figure S7 looking at overall AMR and S6 for specific targets suggest that within the NE there is considerable heterogeneity with much of this small area having a low burden of resistance. Please comment on those factors contributing -and the heterogeneity that is apparent in the data.\n6. The authors substitution of not-detected as half of the limit of detection should be cited – although this is well supported in other chemical wastewater based studies. I do fear however, that this may have influenced some of their analyses.\na. As mcr was identified the least (i.e. negative in 38.3% samples) - how much significance are we seeing in the population modelling as a result of their assumption of not detected being equivalent ot half the limit of detection? What if these were changed to zero? Given its strongest associations with SES factors this would be wise. This is also true for NDM (21.5% non-detects).\n7. The authors do not (to the degree required) – distinguish between association and causation in their discussion. This would be very important particularly as it pertains to the SES factors assessed here – and some caveats through the manuscript would be wise.\na. A multivariate analysis - particularly with such a large dataset to try to disentangle the many observed associations to determine which are independent would be nice.\nMinor issues.\n1. The authors do a good job in the figures – but less so in the text in emphasizing that their data has been normalized for 16S rRNA.\na. They do not also describe the rationle for this step\nb. Please ensure consistency in naming. 16S, 16s rRNA\n2. The authors bounce around between ddPCR and dPCR. Were both used – and if so where and when? Or was only one methodology employed.\n3. There are inconsistencies using proper formatting for genes vs proteins throughout the figures and texts.\n4. I think a spelling mistake here …. “At each split” - or each site? Line 205.\n5. Why are the data not always presented for each of the ARG-targets in each of the figure and analyses? (ie Fig S4 only 7 of the 11 gene targets (many beta-lactamases but not all) - and these differ from Figure S4 where a different 7/11 are presented . Figure S3 has 8/11. Figure S6 11/11.\na. If there is a rational for the presenting only these selective examples please provide that.", "reviewer_3": "The authors measured 11 antibiotic resistance genes by digital droplet PCR in wastewater solids obtained from 163 wastewater treatment plants across the United States. The topic is very interesting, and the work has significant potential for impact. Further clarification on the following aspects is needed for readers to better understand its approach and potential impact.\nData Quality:\n• Line 391: Wastewater concentrations vary significantly both diurnal and weekly. While the diurnal variation can be addressed using composite samples, the impact of weekly fluctuations remains unclear. Limited sampling (2.7 samples per plant) could introduce uncertainties. More details on the sampling and data processing are needed to justify the sufficiency of the average of 2.7 samples per wastewater treatment plant.\n• Line 554: How different are the samples from multiple treatment plants within the same county? This is an important issue to examine before randomly keeping only one single value for these counties to build the RF models.\n• It is not clear why only wastewater solids were analyzed and ARG in dissolved solids or the liquid fraction were excluded in this study. Can we assume the distribution of ARG in the liquid phase and suspended solids to be consistent and stable? Is ARG concentration in the wastewater solids representative of the total concentration in wastewater?\nEffect of airport presence & international travel:\n• How to interpret the effect of airport presence? A sewershed without an airport can either be a rural area or a large/mega city with the main airport just outside of the sewershed boundary. Many major cities have airports situated in neighboring municipalities or on the urban outskirts, rather than within the central urban area. This may explain why so many samples do not feature an airport (Figure 2).\n• Line 271: This conclusion is not supported by the RF results, where the number of airports had a variable importance of zero, for mcr-1 and all other ARGs.\nModel accuracy:\n• Line 214: with R² values ranging from 0.30 to 0.40, can the explanation that a \"significant amount of the underlying variation\" is captured by non-clinical indicators be fully justified? Given the relatively low model performance, more details are needed to assess whether the RF models were trained properly to achieve the best possible accuracy. Key considerations include:\na. The quality of the ground truths: Were the wastewater samples representative of broader trends? More information on sampling strategy would help clarify this.\nb. Input features selection and elimination: How were the 25 input feature selected? Was there a formal feature selection process? How to justify the appropriateness of the selected features?\nc. Line 557: Why were antibiotic usage variables excluded from the RF models? Given their likely impact on ARG concentrations, how reliable are the variable importance results from a model that does not include these key predictors?\nd. Line 204: why were the random forest models tuned using only one parameter, while there are many other parameters that can be tuned? Also, given the complex input-output relationships, how was the decision to use “only 2 variables at each split” justified?\nImportance analysis:\n• A major focus of this study is to understand the factors influencing ARG concentrations. This was explored through preliminary statistical analysis (Figure 2) and variable importance analysis using RF. Given the huge uncertainties with RF, why were alternative statistical/multivariate analysis techniques not considered?"}, "author_response": "We thank the reviewer for their thoughtful and constructive comments. We appreciate the detailed feedback regarding the limitations of our modeling approaches. We have addressed these concerns comprehensively, point-by-point, in our revised manuscript and response to reviewers below, including revising our predictive model to be limited to regions with similar characteristics to our training set, emphasizing the correlational nature of our analysis, and providing additional information on our secondary variables for improved clarity.\n\n[... Author response continues for all points ...]"}, {"round": 2, "reviewer_comments": {"reviewer_2": "I was pleased to re-review the manuscript by Pickering et al., I am pleased to see that it is again being considered for publication – it is an impactful manuscript that is broadly representative of the US population.\nThe authors have done an excellent job of trying to address reviewer’s comments (longer than the paper itself!) in a thorough and thoughtful manner. I am satisfied with these efforts – and the manuscript is both stronger, and easier to read. Most importantly the appropriate caveats and qualifiers are now included.\nA few very, very, very minor comments;\nOccasional areas in which nomenclature for genes is not adhered to (ie vanA in line 2019)\nLines 2610-263 comparing FQ resistant Enterobacterales vs 3rd gen ceph resistant (E. Coli, Kleb, Enterobacter, Citrobacter … these account for >90% of Enterobacterales clinical isolates… this could be rephrased to appear more consistent\nI recommend this for FULL ACCEPTANCE – minor areas can be changed during the page setting step.", "reviewer_3": "I appreciate the authors’ efforts to address the concerns raised during the initial round of review. However, several key issues have not been properly addressed. For example:\n- Sampling limitations: I am unconvinced by the explanation regarding the limited sampling. While I acknowledge the logistical challenges associated with this large-scale study, the core issue lies in the representativeness of the snapshot data. Collecting samples within the same week does not sufficiently address this concern—wastewater characteristics are influenced by weather conditions and local water usage patterns, which can vary significantly across locations, even within the same week. A more robust justification for the representativeness of the snapshot data is required before the detailed results discussion can be considered acceptable. At a minimum, the limitations of the snapshot approach should be explicitly acknowledged and embedded throughout the discussion (though I don’t think this alone would make the approach acceptable.)\n- Random selection for RF modeling: The response regarding the decision to randomly keep a single value for counties with multiple treatment plants in the RF models remains insufficient. This decision needs to be justified by discussing other alternatives.\n- Feature selection: The authors have modified their input feature set (e.g., removing airport presence, adding antibiotic usage variables), but there is no discussion on how these changes affected model performance or sensitivity. This, again, raises concerns about how the selected features are adequately justified.\n- The factors influencing ARG concentrations: The decision to rely entirely on preliminary statistical analysis and RF variable importance for attribution analysis is not well-supported. There are many interpretability methods for models such as XGBoost and neural networks, and the statement that interpretability is lacking in these models is not valid and should be reconsidered. This links back to my earlier concern regarding the reliability of attribution analysis based on RF alone.", "reviewer_4": "The authors present a well constructed response to the initial reviewer feedback. I do not have any additional concerns. This is a very well written article that balances the nuances of WBS with the need to understand drivers in human communities."}, "author_response": "We thank the reviewer for the positive endorsement and another round of careful review. We’ve gone through the manuscript to look for nomenclature mistakes and have rephrased the sentence mentioned for clarity as:\nFor example, in a recent preprint, Goetgeluck et al. cultured bacteria from 12 diverse sewersheds in Atlanta and found that sewersheds with higher proportions of Hispanic, non-Hispanic Asian, non-English speaking residents, and crowded households exhibited higher concentrations of fluoroquinolone-resistant Enterobacterales, and third-generation cephalosporin-resistant bacteria (E. coli, Klebsiella, Enterobacter, or Citrobacter spp.), similar to significant potential determinants found in our study."}, {"round": 3, "reviewer_comments": {"reviewer_5": "I appreciate efforts made both by authors and other reviewers. I have provided my review comments (in red) as a mediating referee in the attachment."}, "author_response": "We appreciate the reviewer’s time and thoughtful consideration of both our manuscript and other reviewers’ comments. Since the mediating reviewer’s suggestions were addressed in our previous revision, this document clarifies where those specific edits and responses can be found in the current text."}]

Decision letter