Wastewater-Based Genomics as a Scalable Framework for Antiviral Resistance Surveillance

1. Introduction

Antiviral resistance is an emerging global health threat that can compromise the long-term effectiveness of therapeutic strategies and undermine pandemic preparedness [1]. While clinical surveillance remains the primary approach for detecting resistance mutations, it is constrained by limited sampling coverage, testing access, and lag times in sequencing and data sharing [2]. These limitations highlight the need for scalable, real-time platforms capable of monitoring antiviral resistance dynamics across entire populations, regardless of healthcare access or testing behavior.

Wastewater surveillance (WWS) has proven to be a powerful tool for community-level tracking of infectious diseases. Throughout the COVID-19 pandemic, WWS was deployed globally to monitor SARS-CoV-2 viral load, detect emerging variants, and identify spike mutations associated with immune escape [3,4,5,6,7,8]. While WWS has been demonstrated for monitoring antibiotic resistance [9,10], the same matrix has yet to be fully leveraged for monitoring resistance to antiviral drugs.

SARS-CoV-2 provides a timely and well-characterized model to evaluate this approach. In response to the pandemic, several antiviral therapies were rapidly developed and deployed. Remdesivir, a viral RNA-dependent RNA polymerase (RdRp) inhibitor, was the first approved antiviral for COVID-19 and received emergency use authorization in May 2020. It was followed by molnupiravir, authorized in late 2021, and nirmatrelvir–ritonavir (Paxlovid), a main protease (Mpro) inhibitor authorized around the same period [11]. In parallel, several monoclonal antibody (mAb) therapies such as bamlanivimab, casirivimab-imdevimab, and sotrovimab were introduced to reduce severe outcomes in high-risk populations [12,13].

Shortly after clinical deployment, resistance-associated mutations were reported in both treated individuals and global genomic datasets. Mutations in the RdRp gene (nsp12), including E802D, C799Y, and M794I, have been associated with decreased susceptibility to remdesivir [14,15,16]. In the Mpro gene (nsp5), mutations such as T169I and A173T emerged under nirmatrelvir pressure and have been linked to resistance phenotypes [14,17]. Similarly, spike mutations such as S:E340D/K/A/V and S:P337L/T reduce the neutralization activity of sotrovimab [12,18,19]. More recently, S:R346T, S:F456L (FLiRT), and S:Q493E (FLuQE) have been linked to broad immune escape, limiting vaccine-elicited neutralization across multiple Omicron sublineages [20,21,22,23].

While these resistance-associated mutations have been characterized clinically, their emergence and dynamics at the population level remain poorly understood, particularly in the absence of widespread and continuous clinical sequencing. WWS has the potential to close this gap by detecting early-phase resistance emergence, regional amplification, and population-scale circulation trends. Here, we used SARS-CoV-2 as a model system to evaluate the utility of wastewater-based genomic surveillance for detecting and characterizing mutations associated with antiviral resistance. We analyzed ~10,000 wastewater samples collected from six public health regions in Ontario, Canada, spanning a 34-month period from October 2021 to July 2024. Our goal was to assess WWS as a scalable, complementary resistance monitoring system, with potential to enhance early-warning capabilities and inform public health responses to antiviral resistance threats.

2. Materials and Methods

Sample collection and processing. As part of the province-wide wastewater surveillance initiative for SARS-CoV-2 variants by the Ontario Ministry of the Environment, Conservation and Parks (MECP) and public health organization and university partners, wastewater samples were collected, preprocessed and sequenced by participating university partners [6,8,24]. Overall, this study includes data from approximately 10, 000 samples collected from across six Ontario health regions over the period of 34 months (October 2021 to July 2024). 24-hr composite samples were collected using different sampling methods including passive torpedo samplers and autosamplers [8,24]. The SARS-CoV-2 nucleic acid quantification, extraction, and tiled-amplicon metagenomic sequencing on Illumina platforms was performed following the protocol used in literature published by the consortium [6,8,24,25].

Data processing. Pair-end sequence data for all samples were quality-filtered and trimmed to remove adapter sequences using cutadapt v3.7 [26]. Thereafter, the remaining high quality reads were mapped to a SARS-CoV-2 reference genome (NC_045512.2) using minimap2 v2.28 [27] that generated an alignment file (bam files) for each of the samples.

Estimating mutation frequency. Mapped genomes were screened with Alcov [28] using the default minimum read depth of 10 for the manually curated experimentally validated mutations that encode resistance to antiviral therapies for SARS-CoV-2 infection described in literature, including those reported in Stanford Coronavirus Resistance Database (CoV-RDB; https://covdb.stanford.edu) [13,14,17,19,29,30,31]. Data were filtered for samples with incomplete metadata. Also, mutations occurring in fewer than 10 samples and frequencies < 1% were filtered.

Statistical analysis. All statistical analyses were performed in RStudio (v1.4.1717). To characterize spatiotemporal patterns of resistance-associated mutations, we applied non-parametric Kruskal-Wallis tests to assess regional differences in frequencies for each mutation. Temporal trends were visualized using LOESS smoothing, with mutation frequencies plotted across the sampling period. Only mutations detected in ≥10 samples at a frequency ≥1% were included in downstream analysis. Descriptive statistics for each qualified mutation included the number of samples detected, median frequency, interquartile range (IQR), region of first detection, peak abundance, and timing of peak signal. p-values <0.05 were considered statistically significant.

3. Results

We assessed 10, 000 SARS-CoV-2 positive wastewater samples for their potential for use in tracking the emergence of resistance of SARS-CoV-2 to antiviral therapies. The frequency of 55 unique mutations that had been experimentally validated to be associated with resistance to antiviral drug therapies in SARS-CoV-2 infection, specifically three main therapies including nirmatrelvir, remdesivir, and sotrovimab were investigated (Supplementary Table S1).

Spatiotemporal characterization of the spike mutation S:E340D associated with sotrovimab resistance. The spike mutation S:E340D, a substitution within the receptor-binding domain (RBD) known to confer resistance to the monoclonal antibody sotrovimab, was detected in 46 wastewater samples collected between October 2021 and July 2024. This mutation was identified in all six Ontario public health regions and met established inclusion criteria (≥10 samples at ≥1% frequency). S:E340D was first detected in the South West region on February 23, 2022, at a frequency of 2.6%, marking one of the earliest wastewater-documented resistance-associated spike mutations in the surveillance network (Table 1). The peak detection occurred in the North East region on February 25, 2022, reaching 5.63%, yet no single region exhibited recurrent elevation over time. Across the dataset, the mutation maintained a median frequency of 1.55% (IQR: 0.85%), with no instances exceeding 10%. Despite broad geographic coverage, Kruskal-Wallis testing revealed no significant regional differences in mutation abundance (p = 0.176), underscoring a pattern of spatially diffuse or low-level circulation (Figure 1). The absence of fixation or persistent elevation suggests that S:E340D emerged transiently under selective pressure, possibly influenced by sotrovimab prescribing practices or shifting population immunity. While no catchment reached dominance by this mutation, its recurrence across multiple sewer sheds and early emergence in the surveillance timeline emphasize the sensitivity of wastewater genomics in capturing early-phase resistance-associated immunoevasive substitutions, even in the absence of clinical sequencing confirmation

Spatiotemporal surveillance of remdesivir resistance associated mutations in wastewater. Abundance of 12 RNA-dependent RNA polymerase (ORF1b) mutations previously implicated in reduced susceptibility to remdesivir were assessed. Of these, five mutations (ORF1b:N189S, S750A, M785I, C790F, and E127V) met predefined inclusion thresholds and were retained for quantitative assessment (Figure 2a–e). The five mutations demonstrated markedly heterogeneous detection patterns in both time and space. ORF1b:S750A was the most frequently detected (n = 142), followed by M785I (135 samples) and C790F (106 samples); N189S and E127V were detected in 21 and 24 samples, respectively. All five mutations were observed in at least four health regions, with S750A, M785I, and C790F exhibiting province-wide distribution. Median frequencies ranged from 1.89% (N189S) to 3.32% (C790F). These values, while subdominant, exceeded detection thresholds consistently and, in some instances, surged. Temporal onset varied across mutations. The earliest detections were observed in the North East region: M785I and C790F emerged in October 2021, both exceeding 5% frequency at first detection (Figure 2c,d, Table 1). In contrast, S750A first appeared shortly after (October 27, 2021) at 1.51% (Figure 2b), while N189S and E127V emerged later, in Toronto (January 2022) and Eastern Ontario (December 2021), respectively (Figure 2a,e). This staggered temporal emergence across geographically distinct catchments suggests either independent introduction events or early cryptic dissemination. Despite modest median frequencies, all five mutations exhibited isolated yet pronounced peaks in specific regions. C790F and M785I reached maximum frequencies of 100% and 99.86%, respectively, in Southwest Ontario, both in 2023, suggesting near-complete dominance of resistance-associated variants within those sewer sheds. N189S peaked at 37.33% in Central East in early 2024, while S750A and E127V reached 7.14% (East) and 9.09% (North East), respectively. These spikes were temporally bounded, and no mutation demonstrated sustained elevation over consecutive surveillance intervals. Statistical testing revealed no significant spatial enrichment for any of the five mutations. Kruskal–Wallis comparisons of mutation frequency across regions yielded p-values ranging from 0.41 to 0.93, suggesting that the observed regional variation reflects stochastic introduction or limited amplification rather than persistent or localized expansion. Notably, even mutations with province-wide detection (e.g., S750A, M785I) failed to show region-specific concentration, indicating diffusion-like patterns of circulation. Taken together, these data demonstrate that remdesivir resistance–associated mutations circulate at low to moderate abundance in Ontario wastewater but appear non-uniformly, transiently, and without consistent spatial anchoring. The detection of such mutations at high frequency in isolated events, approaching or reaching fixation in some catchments could raise important questions about localized treatment pressure, clonal sweeps, or epidemiologically siloed transmission chains. These signals, inaccessible to traditional clinical sequencing frameworks, underscore the unique value of wastewater genomics for detecting resistance emergence below the clinical radar.

Spatiotemporal dynamics of main protease resistance–associated mutations detected in wastewater. Thirty-two mutations described previously to be associated with nirmatrelvir resistance were screened, but only three ORF1a mutations (T3432I, Q3452K, and R3451S) were detected above interpretive thresholds in wastewater collected from six Ontario public health regions (Figure 3a–c). These mutations exhibited distinct geographic and temporal trajectories, including multiple independent episodes of high frequency, and in two cases, complete fixation in wastewater viral populations. Q3452K was the most frequently observed mutation, detected in 264 samples and present in all six health regions. It emerged early, first appearing in North East Ontario on January 10, 2021, at 16.33% frequency, and later reached 100% abundance in the Toronto area on June 21, 2023, representing a transient but complete replacement of circulating lineages in that catchment (Table 1s. Its median frequency was 3.85% (IQR: 6.1%), and Kruskal-Wallis testing confirmed significant regional variation (p = 0.00029), reflecting its concentration in high-population urban zones (Figure 3b). T3432I, though less frequent (38 samples), displayed a similarly punctuated trajectory. It was first detected in Central East Ontario in December 2021 with a strikingly high initial abundance (16.36%) and peaked at 100% in the East region on November 10, 2022. Despite these spikes, its overall abundance was more variable (median = 3.43%, IQR: 6.24%), and regional differences approached but did not reach statistical significance (p = 0.064) (Figure 3a). R3451S was detected in 60 samples across five regions and exhibited a more subdued pattern. First identified in Central West Ontario on December 13, 2021, at 1.5%, it peaked at 100% in South West Ontario in April 2024, suggesting delayed but episodic amplification. However, its lower median abundance (2.74%) and nonsignificant spatial variation (p = 0.263) indicate limited regional fixation or expansion (Figure 3c). Notably, all three mutations exhibited sharp, localized peaks in abundance, frequently exceeding 90%, followed by rapid dissipation. This pattern is consistent with transient selective sweeps or stochastic amplification events within sewer sheds, rather than sustained evolutionary dominance. The early detection of these mutations, often months before peak abundance, demonstrates the sensitivity of wastewater genomics for capturing the leading edge of resistance emergence. Their detection across disparate regions and time points also raises the possibility of convergent resistance evolution under nirmatrelvir pressure in community settings.

Temporal and spatial dynamics of FLiRT and FLuQE mutations in Ontario wastewater. The FLiRT (S:R346T, S:F456L), and FLuQE (S:Q493E) spike mutations previously associated with reduced neutralization by vaccine-elicited antibodies were detected at high prevalence and sustained abundance in wastewater samples collected across six public health regions in Ontario, between October 2021 and July 2024 (Figure 4a–c). S:R346T was the most persistent and dominant of the mutations, detected in 1,539 samples across all six public health regions, representing the highest frequency among all resistance-associated spike mutations in this study. First identified on May 10, 2022, in Toronto area Ontario at 4% frequency (Table 1), it reached 100% in Toronto by September 30, 2022. The mutation sustained a median abundance of 89.12% (IQR: 66.61%) and showed statistically significant regional variation (Kruskal–Wallis p = 0.00001). This mutation exhibited spatial dynamics within the study period (Kruskal–Wallis p = 0.00001) (Figure 4a). This temporal consistency and spatial spread suggest strong selection in populations with widespread vaccine-derived immunity. Similarly, S:F456L was detected in 1,092 samples across all six health regions. It was first observed in South West Ontario on January 11, 2022, at 1.39% frequency, and later peaked at 100% in Central West on August 7, 2023 (median frequency:74.74%; IQR: 48.21%) with significant spatial variability across regions (p = 0.00001). Its detection spanned a 20-month period with consistently elevated signal across multiple catchments (Figure 4b). Prior studies have shown that F456L reduces neutralization by vaccine-elicited sera across multiple Omicron sublineages, including XBB.1.5 and EG.5, and remains resilient in transmission due to minimal fitness trade-offs (Cao et al., 2023; Tamura et al., 2023). Conversely, S:Q493E was detected in 269 samples, with the first occurrence on February 24, 2022, in Toronto, at a low frequency of 1.82%, and peaked on March 08, 2024, in the same region. This mutation had a median frequency of 34.75% (IQR=63.57%) with significant regional differences (p = 0.00285). Although detected less frequently than the FLiRT mutations, S:Q493E exhibited sharp temporal peaks and broad geographic distribution (Figure 4c). The widespread appearance, high frequency, and occurrence of these FLiRT and FLuQE mutations across all monitored regions distinguish them from other resistance-associated mutations in this dataset.

4. Discussion

This study demonstrates the utility of wastewater-based genomic surveillance for detecting and monitoring SARS-CoV-2 mutations associated with reduced susceptibility to antiviral therapeutics and vaccine-elicited neutralizing antibodies. By integrating regionally distributed, longitudinal wastewater sampling with high-resolution amplicon-based sequencing, we captured the presence, temporal patterns, and geographic distribution of functionally important spike and non-spike substitutions at the population level, independent of clinical case reporting or sequencing volume. Among the resistance-associated mutations tracked, those linked to remdesivir and nirmatrelvir targets exhibited divergent patterns. The RdRp mutations associated with remdesivir resistance were observed intermittently and at low-to-moderate abundance, without sustained signal dominance. These patterns are consistent with prior reports indicating that such mutations confer partial resistance but often impose fitness cost in untreated populations [32,33,34]. In contrast, Mpro mutations associated with nirmatrelvir, particularly ORF1a:Q3452K and T3432I, were detected at higher frequencies and displayed regional amplification, and prolonged environmental presence, suggesting local treatment-driven selection. These patterns reinforce the complexity of interpreting antiviral resistance in real-world settings and emphasize the importance of environmental surveillance to detect early ecological signatures of potential drug adaptation.

Notably, remdesivir was introduced in Canada under emergency use authorization in July 2020 [35]. In this study, the earliest wastewater detection of a remdesivir-associated mutation (ORF1b:M785I) occurred in October 2021, shortly after widespread clinical use, while nirmatrelvir-associated mutations (e.g., ORF1a:T3432I) first emerged in late 2021, consistent with early-stage post-authorization exposure [36]. The mutation associated with resistance to monoclonal antibody sotrovimab assessed in this study also emerged in relation to therapeutic and immunization timelines. Sotrovimab received clinical use authorization in late 2021 [37], and its associated resistance mutation, S:E340D, was first detected in wastewater in this study in February 2022, reflecting rapid post-deployment selection. The temporal alignment of mutation emergence with therapeutic rollout supports the utility of wastewater genomics in capturing early signals of treatment-driven resistance evolution.

The FLiRT and FLuQE spike mutations associated with immune escape, (S:R346T, S:F456L, and S:Q493E) showed the most consistent and sustained wastewater signals. These mutations have been documented in experimental studies to impair neutralization by vaccine-elicited antibodies without disrupting ACE2 binding [18,23,38]. S:R346T and S:F456L were detected across all six regions at high frequency and over extended periods, including multiple catchments with near-total prevalence in wastewater viral populations. S:Q493E was less frequently detected but demonstrated episodic dominance and broad geographic distribution. These detection patterns occurred during phases of widespread vaccination [39,40]. The differences in abundance, persistence, and regional concentration between spike mutations and those associated with RdRp. Mpro resistance may reflect differential selection pressures and transmission dynamics. Spike substitutions that reduce vaccine-elicited neutralization can disseminate broadly through population-level transmission, while drug-resistance mutations often arise within treated individuals and may not spread without continued therapeutic exposure [41,42].

While wastewater surveillance does not resolve infection history, treatment status, or clinical outcome, it can reveal the relative scale and spatial patterning of viral mutations over time [8,25]. Importantly, this study highlights the capacity of wastewater genomic surveillance to detect functionally significant substitutions that may not be visible through limited or delayed clinical sequencing. Several of the spike mutations were detected at high abundance in wastewater samples during early phases of regional circulation, underscoring the potential for anticipatory surveillance. Because wastewater integrates signal from all infected individuals, regardless of symptoms, healthcare access, or testing behavior, it enables population-scale monitoring of not only viral evolutionary trends, but also the emergence of antiviral resistance [8,25].

The tracked mutations have documented associations with reduced neutralization or antiviral sensitivity but their contextual relevance may vary depending on co-circulating variants, immunity levels, and therapeutic use patterns not captured in this study [14,32]. Future integration of wastewater genomics with clinical sequencing, prescription data, and serological studies will be necessary to triangulate sources of resistance and evaluate their real-world impact. Nonetheless, this study establishes that wastewater genomics can sensitively and systematically track the emergence and spread of SARS-CoV-2 mutations linked to antiviral resistance and immune evasion. The ability to monitor these mutations in real time across entire populations offers a valuable complement to clinical genomic surveillance, particularly in the context of declining testing volumes and shifting immunity profiles.

5.0. Conclusion

This study demonstrates that wastewater-based genomic surveillance can be used to monitor viral mutations associated with resistance to antiviral therapies. Using SARS-CoV-2 as a model, we tracked the spatiotemporal behavior of mutations previously linked to reduced susceptibility to small-molecule antivirals and neutralizing antibodies. By focusing on substitutions previously linked to reduced susceptibility to remdesivir, nirmatrelvir, and neutralizing antibodies, we show that resistance-associated mutations exhibit distinct patterns of environmental persistence, regional spread, and frequency. This finding supports the application of wastewater genomics for real-time monitoring of antiviral resistance trends, complementing clinical sequencing and enhancing early warning capacity for public health interventions. In addition, it provides a foundation for extending this approach to other viral pathogens requiring ongoing therapeutic monitoring.