Evaluation of a Wastewater-Based Epidemiological Approach to Estimate the Prevalence of SARS-CoV-2 Infections and the Detection of Viral Variants in Disparate Oregon Communities at City and Neighborhood Scales

The prevalence of SARS-CoV-2 in six Oregon communities (some at multiple time points) was estimated using a two-stage cluster sampling scheme as a part of Oregon State University’s Team-based Rapid Assessment of community-level Coronavirus Epidemics (TRACE) project.²⁴ In the first stage, 30 clusters (each comprising one or more census blocks) were randomly selected in each community, with probabilities proportional to the number of housing units (i.e., clusters with a large number of occupied housing units had a higher probability of being selected). Each cluster contained a minimum of 50 occupied housing units. The TRACE project received approval from Oregon State University’s institutional review board to conduct this work.

Starting in 2021 with the last two community sampling events (in Redmond and Corvallis), results from community wastewater testing were used to inform the first stage of sampling. These samples represent a refinement of the TRACE methodology based on evidence that estimates of SARS-CoV-2 concentration in wastewater correlated well at the micro-sewershed level to prevalence within the community. For these two samples, clusters were first grouped into approximately five strata based on micro-sewershed boundaries. Then wastewater data collected $\sim 2–3\mathrm{wk}$ prior to household sampling were used to allocate the 30 sampled clusters to the strata, using optimal allocation to minimize the anticipated standard error of the prevalence estimate. This standard approach in sampling methodology facilitated oversampling from higher prevalence areas in a principled design-based approach and still allowed for approximately unbiased estimates of prevalence to be obtained. With this approach, a total of 30 clusters were still selected; however, a set number of them were required to be from each stratum based on wastewater data. Within each stratum, the clusters selected for door-to-door sampling were still chosen with probability proportional to the number of occupied housing units in the cluster.

In the second stage, field teams used systematic sampling with a random start (i.e., the housing unit at the northwest corner of the cluster) to identify 12 housing units within each of the 30 clusters. The sampling interval (k) used for this stage was specific to each cluster and was calculated by dividing the total number of housing units in the cluster by 12. Field teams were provided with the sampling interval and starting location in advance of data collection. Teams proceeded through their assigned cluster(s) in a serpentine fashion, adhering to the sampling interval. Teams revisited housing units in which no one was home two additional times before the housing unit was replaced. Housing units in which no one agreed to participate were also replaced.

All residents of the selected households, regardless of age, were invited to participate. After obtaining informed consent, field teams interviewed each participant to collect information including name, date of birth, contact information, symptoms, previous positive test result(s), COVID-19 vaccination status, and demographics. Field teams provided test kits to participants and assisted them in self-collecting a nasal swab specimen. Nasal swab specimens were transported to the laboratory for processing and analysis. Test results were mailed to all participants and available to participants via a secure website if they opted for this method. REDCap software was used to securely collect and track data participation and disseminate results of the COVID-19 tests to the participants.²⁵

In total, the TRACE project collected 4,136 self-administered nasal swab samples from residents in 2,521 randomly selected Oregon households in Bend, Corvallis, Eugene, Hermiston, Newport, and Redmond from 30 May 2020 to 14 March 2021. An average of $517\pm 78$ individuals from $315\pm 43$ households ($60\%\pm 14\%$ average household participation rate) participated in each random door-to-door nasal swab sampling event. The prevalence of SARS-CoV-2 within the community at the time of sampling was estimated using one of two approaches based on the number of positive cases detected.

When at least one positive case was detected, a design-weighted approach accounting for the corresponding multistage sampling design (with or without wastewater strata) was used. In this approach, the sampling unit is the housing unit, the sampling weight depends on factors of the design, and the response is the observed positivity among participants in that housing unit. For example, in the original approach without incorporating the wastewater information, the sampling weight of housing unit j in cluster i is shown in Equation 1 where $M_0$ is the total number of housing units in the population; $M_i$ is the total number of housing units in cluster i; $m_i$ is the number of housing units sampled from cluster i (ideally 12); and $h_{ij}$ is the number of participants from housing unit j in cluster i.These weights are used to estimate prevalence using a standard Horvitz–Thompson estimator.²⁶

When zero positive cases were identified through door-to-door sampling in the whole community, a Bayesian approach was taken that combined the observed data with active case counts within the community at the time of sampling. In particular, a beta-binomial model was used where the count of positive individuals was assumed to follow a binomial (n,p) distribution, and the prior for $p$ is assumed to be beta ($\mathrm{\alpha }$, $\mathrm{\beta }$). Values for $\mathrm{\alpha }$ and $\mathrm{\beta }$ were obtained via optimization by setting the fifth quantile of the beta distribution equal to the number of cases reported to the Oregon Health Authority in the previous 2 wk divided by the population size, and the 50th quantile equal to the fifth quantile divided by the proportion of cases thought to be symptomatic based on previous data.

A noninformative (uniform [0,1]) prior was also explored, and the approach was not sensitive to the choice of prior. The posterior distribution for $p$ also takes the form of a beta and can be used to estimate prevalence. Imperfect test sensitivity and specificity were accounted for in both approaches. Prevalence estimates were quickly shared with the local public health officials partnering with the TRACE project.

Weekly COVID-19 case data for the participating communities were obtained from ZIP code reported data from the Oregon Health Authority website.²⁷ For weeks when $<10$ cases occurred in a given ZIP code, the Oregon Health Authority reported the case count as “1–9 cases.” In these instances, to minimize error, a value of 5 cases was substituted in the calculations. The ZIP codes included for each community were as follows: Bend (97701, 97702, 97703, 97707, 97708), Corvallis (97330, 97331, 97333), Eugene (97401, 97402, 97403, 97404, 97405, 97447, 97478), Hermiston (97838), Newport (97365, 97366), and Redmond (97756). Reported test results with location data for the city of Newport were obtained with permission from the Lincoln County Health and Human Services Department, and public health officials anonymized the data by aggregating according to micro-sewershed prior to analysis.

Several wastewater concentration methods were tested to determine which method resulted in the highest viral recovery, as estimated by the maximum mean SARS-CoV-2 signal recovered from replicate composite WWTP samples. The WWTP sample was collected from Newport, Oregon, during the June 2020 sampling event. Two concentration methods, electronegative membrane filtration (EMF) and centrifugal ultrafiltration (CU) were evaluated. In addition, four iterations of EMF with various amendments were tested (EMF Methods A–D). All methods were conducted in triplicate.

For the EMF method, $30\mathrm{mL}$ of WWTP influent composite sample was mixed gently on a stir plate with a magnetic stir bar. Each of the EMF method iterations differed in the second step, when amendments were made to the sample. In EMF Method A, 0.01 N hydrochloric acid (HCl) was added to titrate the pH below 3.5 and magnesium chloride (${\text{MgCl}}_2$) was added to make a final concentration of $25\text{\hspace{0.17em}mM}$. In EMF Method B, 0.01 N HCl was used to titrate the pH below 3.5 but no ${\text{MgCl}}_2$ was added. In EMF Method C, ${\text{MgCl}}_2$ was added for a final concentration of $25\text{\hspace{0.17em}mM}$, but the sample was not acidified. In EMF Method D, no additions were made.

All samples were gently mixed for an additional 5 min. After mixing, the entire sample was vacuum filtered through a $0.45\text{-}\mathrm{\mu }\mathrm{m}$ pore size, $47\text{-mm}$ diameter, mixed cellulose ester electronegative membrane (Catalog no. 7141-104; Whatman). The filter was stabilized in $1\mathrm{mL}$ of DNA/RNA Shield (Zymo Research, Irvine, CA) and stored at $–20º\mathrm{C}$ until RNA extraction and reverse transcriptase droplet digital PCR (RT-ddPCR) could occur.

In the CU method, $40\mathrm{mL}$ of wastewater was added to a Centricon Plus-70 centrifugal ultrafiltration device (Catalog no. UFC701008; Merck Millipore) with a 10-kDa cutoff. The ultrafiltration device was centrifuged for 20 min at $\mathrm{3,220}\times g$. The filtrate was discarded, and the concentrate spun out of the ultrafilter back into the concentrate cup at $\mathrm{1,000}\times g$ for 2 min. The concentrate was removed and placed in $1\mathrm{mL}$ of DNA/RNA Shield.

All WWTP influent samples comprised 24-h time-weighted composites taken prior to primary treatment. The characteristics of each WWTP are given in Table 1. Samples were collected at the time of each random door-to-door nasal swab sampling event, as well as an additional one to three times per week for 6–11 months from April 2020–May 2021.

Twenty-two pump stations serving Newport, Oregon’s Vance Avery WWTP sewershed, were sampled hourly for 24 h from 20–21 June and 11–12 July 2020 (Figure 1). Some pump stations were sampled twice within a single weekend, for a total of 52 pump station samples. The total population of the Newport sewershed was 10,853, with the populations served by each micro-sewershed ranging from 15 to 9,426 persons. The characteristics of each micro-sewershed and the distribution of samples are given in Table 2.

For both WWTPs and pump stations, the 24-h composites consisted of hourly samples and were kept on ice during sampling. Samples collected on or before 31 July 2020 were frozen in $200\text{-mL}$ aliquots and stored at $-80º\mathrm{C}$ for up to 33 d prior to concentration (median 6 d). After thawing in cold water baths, samples were concentrated using electronegative filtration, as previously described.²⁸ Briefly, all samples collected on or before 10 July 2020 were acidified to a final pH of 3.5 and ${\text{MgCl}}_2$ was added to a final concentration of $25\text{\hspace{0.17em}mM}$. Samples ($30–40\mathrm{mL}$) were vacuum filtered through a $0.45\text{-}\mathrm{\mu }\mathrm{m}$ pore size, $47\text{-mm}$ diameter mixed cellulose ester electronegative filter (HAWP; Millipore).

Influent samples collected after 10 July 2020 were filtered with no amendments given that preliminary data showed improved viral recovery with unamended wastewater (Figure S1). In addition, influent samples collected after 31 July 2020 were neither frozen nor amended prior to filtration. In these samples, filtration occurred within 8 h of sample collection.

Once filtration was complete, the electronegative membranes were placed into $2\text{-mL}$ tubes containing either $0.7\text{-mm}$ garnet or $0.5\text{-mm}$ glass beads, stabilized in $1\mathrm{mL}$ of DNA/RNA Shield (Zymo Research), and frozen until analysis. Preliminary experiments showed no difference between these homogenization methods (Figure S2). Field blanks of deionized water were processed with every batch of samples.

Participant nasal swab samples were analyzed at the Oregon Veterinary Diagnostic Laboratory using the TaqPath COVID-19 Combo Kit (Applied Biosystems), in accordance with the manufacturer’s instructions for use as required by the Emergency Use Authorization under strict biosafety level 2 (BSL2) conditions. Nucleic acid isolation was performed using the MagMax Viral/Pathogen II Nucleic Acid Isolation Kit (Applied Biosystems/ThermoFisher). Briefly, $200\mathrm{\mu }\mathrm{L}$ of transport medium from the swab sample was added to a single well of a KingFisher Deepwell 96-well plate containing $5\mathrm{\mu }\mathrm{L}$ of Proteinase K. Each 96-well plate held 94 participant samples and 1 negative control well containing water. The last well was left empty to allow for the positive control to be added during the (RT)-PCR detection step. After sample addition, nucleic acid magnetic beads were resuspended, $10\mathrm{\mu }\mathrm{L}$ was added to $265\mathrm{\mu }\mathrm{L}$ of binding solution and then added to the wells. MS2 phage control ($5\mathrm{\mu }\mathrm{L}$) was added to all wells as an extraction control and processed on a KingFisher Flex magnetic particle processor. Purified nucleic acids were eluted in $50\mathrm{\mu }\mathrm{L}$ of MagMax elution solution. Eluted nucleic acid was stored at $-80º\mathrm{C}$ unless RT-PCR was run within 2 h of extraction.

Detection of SARS-CoV-2 viral RNA was performed using the TaqPath RT-PCR COVID-19 Kit on a 7500 Fast Real-Time PCR Instrument (Applied Biosystems). Reactions were run in multiplex with primers and probes specific for three gene sequences specific to SARS-CoV-2: ORF1ab, N Protein, and S Protein. A primer and probe set was included to detect the MS2 phage added during the initial sample processing as an internal control to verify RNA extraction. The thermal protocol for the one-step real-time PCR was as follows: 2 min at 25°C uracil $N$-glycosylase incubation, 10 min at 53°C reverse transcription, 2 min at 95°C activation, followed by 40 cycles of 3 s at 95°C denaturation and 30 s at 60°C anneal/extension/detection. Results were analyzed using SDS Software (version 1.5.1) and interpreted using COVID-19 Interpretive Software (version 1.2; Applied Biosystems).

Wastewater samples were homogenized with either $0.7\text{-mm}$ garnet beads or $0.5\text{-mm}$ glass in DNA/RNA Shield using either a Qiagen TissueLyser (Qiagen Inc.) or a BioSpec Mini-Beadbeater 16 (BioSpec Products, Inc.) for 2 min. Beads and debris were pelleted by centrifugation at 12,000 rcf for 1 min. The lysate was transferred from each tube to a 96-well plate, and $200–400\mathrm{\mu }\mathrm{L}$ was extracted using the MagMAX Viral/Pathogen kit on a KingFisher Flex automated instrument (ThermoFisher Scientific), as described above. Purified RNA was eluted in $50\mathrm{\mu }\mathrm{L}$ of MagMAX Elution Solution. Extraction recovery from RNA extraction step was quantified with a commercial standard (Exact Diagnostics), and extraction blanks were included with every run. RNA was stored at $-80º\mathrm{C}$ until analysis.

Two SARS-CoV-2 RNA targets (N1/N2) and an internal control (Human RNase P) were measured via RT-ddPCR using a commercial triplex assay (2019-nCoV CDC ddPCR Triplex Probe Assay; Catalog no. 12008202; Bio-Rad) using the One-Step RT-ddPCR Advanced Kit for Probes and run on a QX-200 ddPCR system (Bio-Rad). The primer and probe sequences were published previously.²⁹ RT-ddPCR was chosen over RT-qPCR owing to its superior sensitivity,^30,31 robustness against inhibitors,³² and its wide application in WBE methodologies.^17,33,34

Reactions were partitioned into droplets using an automated droplet generator (ADG). Twenty-two microliter reactions were prepared with $5.5\text{-}\mathrm{\mu }\mathrm{L}$ template RNA, whereas the ADG partitioned only $20\mathrm{\mu }\mathrm{L}$, yielding an effective template volume of $5\mathrm{\mu }\mathrm{L}$. Each reaction had an average of 12,657 droplets [$\text{standard deviation\hspace{0.17em}}\left(\text{SD}\right)=\mathrm{1,783}$]. Commercially prepared RNA standards and negative controls were included on each extraction plate and ddPCR plate (cat. no. COV019 and COV000; Exact Diagnostics). All samples and controls were analyzed in duplicate.

The one-step thermal cycling conditions were as follows: reverse transcription at 50°C for 60 min; enzyme activation at 95°C for 10 min; 40 cycles of denaturation at 94°C for 30 s followed by annealing/extension at 55°C for 60 s; enzyme inactivation at 98°C for 10 min; and, last, a 4°C hold for droplet stabilization, for a minimum of 30 min to a maximum of overnight. Finally, the amplification in the droplets was determined using the Bio-Rad droplet reader. All assay conditions were performed as specified in the Bio-Rad assay protocol.³⁵

The quality controls for the wastewater RT-ddPCR method included field blanks, extraction blanks, negative control reactions [containing only human genomic DNA (gDNA)], positive control reactions (containing synthetic RNA of SARS-CoV-2 assay targets and human gDNA), and no-template controls (NTCs). The results of all quality control (QC) reactions are summarized in Table S1. Only positive detections which passed quality assurance/QC were included in this study.

Wells with $<\mathrm{6,000}$ droplets were omitted. Sample data with positive reactions were accepted only if the corresponding extraction blank and field blanks, as well as the PCR negative and no-template controls, were all negative for the N1/N2 targets. When averaging sample data across replicates, a value of one-half the sample-specific limit of detection was substituted for nondetects. Reactions were regarded as positive if three or more droplets per well amplified in either target. Droplet clusters were manually called for each target using the QuantaSoft Analysis Pro software (version 1.4; Bio-Rad). All other analyses were conducted in R (version 4.1.0; R Development Core Team) with Rstudio Desktop (version 1.4.1717).³⁶ Spatial graphics were created using the sf,³⁷ ggplot2,³⁸ gtable,³⁹ and ggmap⁴⁰ packages.

The N1 and N2 markers exhibited generally good agreement in the wastewater samples. The markers were concordant in 91.4% of reactions: 4.1% of reactions ($n=954$) had a positive detection in N1 only (using a threshold of three positive droplets per reaction), and 4.5% of reactions were positive in N2 only. Quantitatively, the markers were also reasonably well aligned according to a simple linear model [$\text{slope}=0.94$, with N2 as the response variable, adjusted $r^2=0.97$, $\text{root mean square error\hspace{0.17em}}\left(\text{RMSE}\right)=0.58$ copies per reaction]. Accordingly, the N1 and N2 data were averaged together for all reported concentration values.

For limit of blank (LOB) determination, 104 reactions were run with the triplex assay across three plates using the Exact Diagnostics Negative Control as a template. Due to the nonnormal distribution of the results, the LOB was determined using a nonparametric (rank order) method with a false-positive probability ($\mathrm{\alpha }$) of 0.05.

The limit of detection (LOD) was predicted for each target (N1 and N2) using Equation 2, as follows:where $stDev\left(CopiesPerReactio{n}_{LOB}\right)$ is the standard deviation (Dev) of the copies per reaction from the LOB assay. The predicted LOD was subsequently tested by running 60 test reactions at concentrations of 4 and 12 copies per reaction of target using the Exact Diagnostics Standard for SARS-CoV-2.

To determine the viral recovery efficiency of the wastewater processing method used in this study, a spike-in experiment was performed using bovine coronavirus (BCoV) as a surrogate for SARS-CoV-2. BCoV solution was prepared from freeze-dried Calf Guard cattle vaccine (Zoetis). After rehydrating in 3 mL of sterile diluent provided by the manufacturer, the BCoV solution was divided into $100\text{-}\mathrm{\mu }\mathrm{L}$ aliquots and stored at $–20º\mathrm{C}$. To use, the BCoV stock solution aliquot was thawed on ice and vortexed thoroughly; each aliquot was used for a maximum of two freeze–thaw cycles. BCoV stock solution was spiked into wastewater samples at a ratio of 1/1,000 (vol/vol) before the concentration step.

To determine the concentration of the BCoV stock solution, $10\mathrm{\mu }\mathrm{L}$ of the stock was spiked into $390\mathrm{\mu }\mathrm{L}$ of phosphate-buffered saline (PBS), and $200\mathrm{\mu }\mathrm{L}$ of that mixture was extracted following the same protocols used for wastewater sample extractions as described in the main text. One extraction blank, prepared with PBS, was included on each plate as an RNA extraction contamination control. The extracted RNA was then serially diluted (1:10) in nuclease-free water for six dilutions and ran in duplicate using a previously established BCoV assay, following the one-step RT-ddPCR procedure.⁴¹ Stock concentration of BCoV was around 230,000 gene copies $\left(\text{gc}\right)/\mathrm{\mu }\mathrm{L}$.

Process efficiency (i.e., viral recovery) was calculated by dividing the final quantity of BCoV measured in wastewater samples by the quantity of BCoV spiked to each wastewater sample before concentration (Equation 3). Nonspiked wastewater samples were also quantified for BCoV to assess background concentrations.

For cDNA synthesis, $11\mathrm{\mu }\mathrm{L}$ of lysate from positive wastewater or participant samples, were used for single-strand cDNA synthesis, using the Thermo Superscript IV kit with the following modifications: No host gDNA/RNA removal steps were performed, no RNase H step was performed, and the reverse transcriptase incubation step at 50°C was increased from 10 min to 30 min.

The cDNA was used for amplification and sequencing, using the Swift Amplicon SARS-CoV-2 Panel (AL-COV48) together with Swift Amplicon Combinatorial Dual Indexed Adapters (AL-S1A96, AL-S1BA96, AL-S2A96, AL-S2B96). The Swift Biosciences protocol was followed, except that, after optimization experiments, the volume of the G1 reagent was reduced to 25% of recommended level. This resulted in an increase in read coverage for wastewater samples of 2- to 8-fold for the experiments described in this study (with the release of version 2 of the Swift primer set, the volume of this reagent was no longer reduced). Individual libraries (30–96 samples) were quantified by fluorescence, normalized, and pooled for $2\times 150\text{\hspace{0.17em}bp}$ sequencing on a lane of an Illumina HiSeq3000 sequencer.

Except for some initial sequencing experiments, libraries from individuals were prepared on different days and sequenced in different lanes than samples from wastewater to reduce the possibility of contamination of sequences from low titer wastewater samples with sequences from high titer individual samples. Pooled libraries from wastewater samples were often run on two to three lanes of the Hi-Seq 3,000 to increase read depth.

Wastewater sequencing on the HiSeq3000 produced a median percentage of reads mapped of 0.86% compared with a median of 21.7% for nasal swab samples. In general, wastewater samples with an RNA concentration of ${\log }_{10}>4.0$ gene copies of N1/N2 per liter of wastewater (gc/L) reliably produced usable amounts of sequence. Samples with lower RNA concentrations produced variable amounts of sequence, and those ${<\log }_{10}3.0\text{\hspace{0.17em}gc}/\mathrm{L}$ were routinely unsuccessful.

After demultiplexing, Illumina primer sequences were trimmed using BBDuk (v.38.96), sequences were aligned to the reference sequence (NC_045512.2; Wuhan-Hu-1) using BWA-MEM (v.0.7.17-r1188), and the SARS-Cov-2 primers were removed using the Swift Biosciences Primerclip package. Genome Analysis Toolkit software (GATK) (v.4.2.0.0) was used to identify variants compared with the reference sequence and to count the numbers of reference and variants reads at each single nucleotide polymorphism (SNP) site. The ploidy was set to 4 and the down-sampling function was not employed. By not using the down-sampling function, all sequence reads that passed standard QC parameters were used in making variant calls. When required, an Integrated Genomics Viewer (IGV) was used to manually inspect sequence alignments and variant calls.

To identify viral genotypes represented in wastewater RNA by multi-locus sequencing typing, sequences from individual samples from Oregon were used to define multi-locus genotypes as a set of polymorphic sites that were unique to each variant and were not shared with any other variant known at the time (Table S2).

To estimate the fraction of each variant present in the wastewater RNA sequences, the number of reads supporting each variant were summed across all variant-specific SNP sites and divided by the total number of reads spanning those SNP sites. At SNP sites where the total number of reads was $>100$, both the variant read number and total read number were scaled down to a total read number of 100 prior to the estimation calculation. At SNP sites with $>100$ total reads, the actual number of reads was used, thus decreasing the weight of those counts proportionally to the read coverage.

After this estimation was conducted for each variant, the fractions were summed for all identified variants. If the sum was $>0.7$, it was assumed that the identified variants comprised all of the RNA molecules present and each fraction was divided by the total to normalize the total to 1.0. If the initial sum was $<0.7$, it was assumed that the difference was comprised of RNA molecules from unidentified variants, and the fraction of RNA attributable to the unidentified variants was set to the difference between the total and 1.0. If the sum was $>1:5$, the SNP data were manually reviewed to identify and remove artifacts (miscalled 1-bp indels) and remove common SNPs that were shared between the variants present to prevent double-counting. The custom software to conduct these calculations was implemented in R and packaged into a SnakeMake pipeline.⁴² All individual sequences were deposited in the data sharing platform Global Initiative on Sharing All Influenza Data (GISAID; see Table S3 for accession numbers) and all wastewater sequences were deposited in NCBI’s short read archive, under BioProject PRJNA719837.^43,44

To compare the accuracy of SARS-CoV-2 wastewater concentrations and clinically reported COVID-19 cases as estimators of community prevalence, Monte Carlo simulations were performed to account for the uncertainty in the point estimates for each sampling event. Clinically reported COVID-19 cases included cases identified through the TRACE random door-to-door sampling events as well as cases reported through standard health surveillance efforts by the Oregon Health Authority. The Monte Carlo simulations accounted for the uncertainty in the point estimates for each sampling event with the underlying assumption that the community prevalence based off the TRACE random door-to-door sampling events were the closest to the ground truth.

For each simulation, a new wastewater concentration (${\log }_{10}\text{\hspace{0.17em}gc}/\mathrm{L}$) was redrawn for each community sample from a Gaussian distribution with mean equal to the point estimate and standard deviation equal to the standard error of the point estimate. Similarly, a new prevalence was redrawn for each community from a method-specific distribution: a truncated Gaussian distribution with mean equal to the point estimate, standard deviation equal to the standard error of the point estimate, and lower truncation bound equal to zero when the design-based estimator was used; or a Beta distribution with shape and scale parameters estimated using optimization to fit the posterior mean and 95% credible interval when the Bayesian estimator was used.⁴⁵ The reported COVID-19 case numbers were not redrawn because these were the publicly available standard. For each simulation, a simple linear regression model was fit using the new wastewater concentration draw to estimate the ${\log }_{10}$ of the new community prevalence draw, and a separate model was fit using the ${\text{Log}}_{10}$ of the observed COVID-19 case count to estimate the ${\log }_{10}$ of the new community prevalence draw. For each model, the slope, intercept, and $R^2$ value were recorded. In total, 10,000 simulations were performed.

The LOB for N1 and N2 was 2.0 and 4.2 copies per reaction, respectively. Both of these values were below the three-droplet threshold; only reactions with one or two droplets yielded copy numbers at that level, which provided confidence in the chosen threshold for calling positive reactions. All LOB reactions were below the positive threshold for N1, and only 4/104 nontarget reactions had three or more droplets in N2, which was a false-positive rate of 4%.

The predicted LOD based on the LOB results were 4 and 12 copies per reaction for N1 and N2, respectively. For an LOD estimate to be valid, $>95\%$ of test reactions at the predicted LOD value need to amplify above the LOB. The LOD of N2 was confirmed to be 12 copies per reaction, as 58/60 (97%) test reactions at that concentration had copy numbers above the N2 LOB. The N1 LOD was somewhere between 4 and 12 copies per reaction: at 12 copies per reaction, all 60 reactions amplified above the N1 LOB, but at 4 copies per reaction, 13/60 (22%) reactions amplified below the N1 LOB. Using a parametric method (which is an imperfect estimate in this case because the test reaction data were not normally distributed), the N1 LOD was estimated to be 8 copies per reaction. The process efficiency based on BCoV recoveries was 57$\pm 4\%$. BCoV was not detected in nonspiked wastewater samples.

Over the course of 10 months, from May 2020 through March 2021, eight sampling events were conducted in six Oregon communities by the TRACE project to determine the COVID-19 prevalence within those communities. The selected communities represented a diverse cross-section of Oregon, ranging from a small coastal commercial fishing community (Newport), to mid-sized and large university communities in the temperate Willamette Valley (Corvallis and Eugene), mid-sized and large arid high desert communities (Redmond and Bend), and a small agricultural community in eastern Oregon (Hermiston).

Community prevalence results are summarized in Table 3. The response rates ranged from 38% to 71%, with an average response rate of $60\pm 14\%$. The estimated SARS-CoV-2 prevalence, including both symptomatic and asymptomatic infections, ranged from 8/10,000 to 1,687/10,000 persons. The SARS-CoV-2 concentrations in the influent wastewater corresponding to the door-to-door nasal swab samples ranged from 2.92 to $5.13{\log }_{10}$ gene copies per liter of wastewater (gc/L).

The majority of the prevalence sampling events took place in the absence of significant precipitation (Table S4). The only exception was in Newport on 20 June 2020. During this 24-h period, 3.56 mm of rain fell, which increased the wastewater flow rate by 5% compared with 11 July 2020, a 24-h period in which no rain fell.

The SARS-CoV-2 concentrations in the influent wastewater of these same cities over a 6- to 11-month period ranged from nondetect to $5.58{\log }_{10}\text{\hspace{0.17em}gc}/\mathrm{L}$ (Figure S3). The wastewater SARS-CoV-2 concentrations showed a moderate positive correlation (Pearson’s $r=0.71$) when compared with the log to the base 10 of reported cases per 10,000 persons (to normalize for differences in population size). The accuracy of estimating reported COVID-19 cases using wastewater SARS-CoV-2 concentrations was also moderate [$\text{root mean square logarithmic error\hspace{0.17em}}\left(\text{RMSLE}\right)=0.14$ and mean absolute percentage error $\left(\text{MAPE}\right)=0.29$], with wastewater concentrations differing by up to 62-fold representing the same number of reported COVID-19 cases per 10,000 persons (Figure 2A; Tables S6 and S8). Similar correlation strengths and levels of accuracy were observed with each city when analyzed individually (Figure S4).

When the wastewater SARS-CoV-2 concentrations were compared with the TRACE-informed prevalence estimates that include both symptomatic and asymptomatic cases (Figure 2B, Table 3), the positive correlation was much stronger (Pearson’s $r=0.96$), and the accuracy was much higher, as demonstrated by the lower RMSLE (0.03) and MAPE (0.17) values. This correlation also suggests that the detection limit of our WBE method was 3 infections/10,000 persons. Compared with wastewater SARS-CoV-2 concentrations, reported COVID-19 cases were more weakly correlated with estimated prevalence (Pearson’s $r=0.85$) and had a lower accuracy ($\text{RMSLE}=0.05$ and $\text{MAPE}=0.31$) than wastewater SARS-CoV-2 concentrations (Figure 2C; Tables S6 and S8).

In addition, Monte Carlo simulations demonstrated that wastewater SARS-CoV-2 concentrations were significantly better than reported COVID-19 cases counts at estimating COVID-19 prevalence, even after accounting for uncertainty inherent in the wastewater and prevalence estimates (Figure 3; Tables S6 and S8). The median $R^2$ for the wastewater SARS-CoV-2 concentration model was 0.82 compared with 0.71 for the reported-cases model, based on the Monte Carlo simulations. A Wilcoxon rank-sum test, which determines whether the two distributions for $R^2$ have a location shift of 0 vs. the alternative that the wastewater SARS-CoV-2 concentration $R^2$ distribution has a positive shift, gave a $p$-value of $<0.0001$. Thus, even after accounting for the uncertainty inherent in the wastewater SARS-CoV-2 concentration and community prevalence estimates, the difference between the wastewater SARS-CoV-2 concentration and reported-cases median $R^2$ values was significant, and the wastewater SARS-CoV-2 concentration had a larger median $R^2$.

For each simulation, the two nonnested regression models were directly compared using Vuong’s test, which is a likelihood-ratio–based test for model selection using the Kullback–Leibler information criterion.^46,47 Vuong’s test revealed that the wastewater SARS-CoV-2 concentration model fit significantly better than the reported-cases model ($p<0.05$) 53% of the time. The reported-cases model fit significantly better only 2% of the time; the models could not be distinguished on the remaining occasions. Although Vuong’s test cannot determine whether the preferred model is the true model, the high $R^2$ value (0.91 using the observed data) and low RMSLE (0.03) and MAPE (0.17) are indicative of a good fit.

To identify COVID-19 “hotspots” at the neighborhood-scale in Newport, the wastewater SARS-CoV-2 concentrations were quantified from wastewater samples collected at 22 pump stations located throughout the small community of 10,853 persons (Figure 1). These pump stations divided the community into neighborhood-scale micro-sewersheds with estimated populations ranging from 15 persons at the 10th St. pump station up to 9,426 persons at the Northside pump station (Table 2). These wastewater samples were collected during the same time period as the two Newport random nasal swab prevalence sampling events.

During the peak of a Newport outbreak in mid-June, 2020, the random nasal swab sampling event estimated a community COVID-19 prevalence of 3.4% (Table 3). During this time, SARS-CoV-2 was detected in wastewater samples from 13/22 micro-sewersheds, with concentrations ranging from $3.16\text{\hspace{0.17em}to\hspace{0.17em}}5.25{\log }_{10}\text{\hspace{0.17em}gc}/\mathrm{L}$ (Figure 4; Table S5). Three weeks later, in July 2020, the estimated COVID-19 prevalence in Newport had decreased to 0.6%. The number of micro-sewersheds with detectable levels of SARS-CoV-2 in their wastewater also decreased to 5, with concentrations ranging from $3.36\text{\hspace{0.17em}to\hspace{0.17em}}4.71{\log }_{10}\text{\hspace{0.17em}gc}/\mathrm{L}$. In addition to decreased micro-sewershed wastewater SARS-CoV-2 concentrations, the decrease in community prevalence also corresponded with marked decreases in nasal swab positivity rates and clinically reported COVID-19 cases (Figure 4; Tables S5 and S7).

Wastewater SARS-CoV-2 concentrations collected from Newport’s 22 micro-sewersheds were significantly correlated with the positivity rate of random swab samples [Pearson’s product-moment $\text{correlation}=0.64$; 95% confidence interval (CI): 0.32, 0.83, $\text{df}=22$; $p<0.001$] and reported cases (Pearson’s product-moment $\text{correlation}=0.61$; 95% CI: 0.39, 0.76; $\text{df}=44$; $p<0.001$). Similarly, the presence of SARS-CoV-2 in wastewater was significantly correlated with reported cases within each micro-sewershed [Fisher’s exact test; $\text{odds ratio\hspace{0.17em}}\left(\text{OR}\right)=24.3$; 95% CI: 3.1, 341; $p<0.001$, $n=33$]. In addition, the presence/absence of SARS-CoV-2 was also significantly correlated between wastewater and random nasal swabs within each micro-sewershed (Fisher’s exact test; OR could not be calculated; 95% CI: 2.16, $\infty$; $p=0.003$; $n=22$).

Using Multi-Locus Sequence Typing,⁴⁸ two distinct viral variants, a B.1.399 lineage variant designated NA and a B.1 lineage variant designated NB, were detected both in the wastewater and among the individuals who tested positive during the June 2020 Newport sampling event (Table S2). The B.1.399/NA consensus most closely matched sequences found in Europe and then in California in March and April 2020, respectively. The B.1/NB consensus sequence most closely matched that from Yakima County, Washington, on 29 April 2020, and sequences found in Europe in March 2020.

During the June sampling event (from 569 samples), B.1.399/NA was detected in 77% of the individuals who tested positive (10/13), whereas B.1/NB was detected in 23% of the individuals who tested positive (3/13) (Figure 5, Table 4). In the July sampling event, two individuals tested positive (from 550 samples). Although both samples both yielded low coverage sequence data, they were identified as B.1.1.291 by the international database, GISAID and, thus, were different from the individuals from the June sampling event.

In the wastewater samples, the fraction of viral RNA accounted for by each variant was estimated from the total fraction of variant reads summed across all SNP positions specific to that variant. During the June sampling event, the viral variant distribution in the WWTP influent was dominated by B.1.399/NA, accounting for 70% of the reads, whereas B.1/NB accounted for a minority of the reads at 4% (Figure 5, Table 4). This mirrors the observations made through the random nasal swab assays (77% B.1.399/NA and 23% B.1/NB).

B.1.399/NA was also dominant in 11/12 positive micro-sewershed wastewater samples collected. During the June sampling event, B.1.399/NA accounted for 30–98% of the total viral sequence reads across the micro-sewersheds, whereas during the July sampling event it accounted for 48–92% of the viral sequence reads across the micro-sewersheds (Table 4). In addition, although the abundance of B.1/NB was always in the minority, ranging from 0% to 52% of total virus sequence reads during the June sampling event and 0% to 11% of total virus sequence reads during the July sampling event, it was always detected in the wastewater of micro-sewersheds, when B.1/NB was detected among individuals via random nasal swabs (Figure 5, Table 4).The only exception to this relationship was the Bayfront micro-sewershed. In Bayfront, B.1/NB was not detected in the wastewater, and the single positive individual discovered in Bayfront carried B.1/NB.

Finally, from the weeks of 8 June to 30 November 2020, five SARS-CoV-2 variants were detected consistently in wastewater samples from the Newport WWTP (Figure 6A), representing at least 5% of the sequence reads in samples from at least 3 wk. Variant B.1.399/NA was the dominant variant during the weeks of 8 June to 6 July, with lesser amounts of B.1/NB detected. During the weeks of 20 July and 27 July, a B.1.369 subvariant, designated EE (Table S2B), was most abundant, and it was also detected during the weeks of 29 June and 3 August. From 3 August onward, a B.1.2 subvariant, designated FF, was the most abundant variant detected. The detection of these variants in Newport broadly mirrored trends across the state of Oregon observed in individuals whose positive samples were sequenced and deposited in GISAID (Figure 6B; Figure S5).

WBE has been widely used to quantify wastewater SARS-CoV-2 concentrations during the COVID-19 pandemic to provide a relative sense of viral burden in a community and track its trend over time.^2,49,50 However, the moderate correlation between wastewater SARS-CoV-2 concentrations and reported COVID-19 cases, similar to those observed in this study, has resulted in the perception that WBE is limited in its ability to estimate community COVID-19 infections.⁵¹ It has been hypothesized that this moderate correlation can be attributed to a variety of factors related to the wastewater sample, including virus decay in sewage conveyance systems, variability in sampling and concentrating techniques, and variability in the magnitude and duration of shedding by infected individuals.^6,11,16 It should be noted that rainfall events were largely absent during the wastewater collection periods of this study. Thus, dilution of wastewater virus concentrations by rainfall events was not an important source of the observed variability in wastewater SARS-CoV-2 concentrations.

Although the correlation between wastewater SARS-CoV-2 concentrations and COVID-19 cases (as reported by health departments) has been moderate,^51,52 this study demonstrated a very strong correlation between wastewater SARS-CoV-2 concentrations and estimated COVID-19 prevalence. With an average response rate of $60\pm 14\%$, the response rate of this study was in the upper end of the range reported by other studies doing a door-to-door health survey (20–50%), lending strength to the correlation analyses.^53–55 In addition, Monte Carlo simulations indicated that wastewater SARS-CoV-2 concentrations were significantly better ($p<0.05$) than reported COVID-19 cases at estimating COVID-19 prevalence in a community. The COVID-19 prevalence estimates minimized many of the uncertainties inherent to reported COVID-19 cases, including asymptomatic individuals and limited access to testing. This suggested that the weaker correlations between wastewater SARS-CoV-2 concentrations and reported COVID-19 cases had less to do with the inherent variability within the wastewater measurements and more to do with the inherent variability in the reporting of COVID-19 cases. In addition, the inherent variability in reported COVID-19 cases is expected to become even more complex with the anticipated unreliability of self-reporting positive results from at-home testing kits.

Thus, when compared with a metric with much smaller uncertainties, such as the prevalence estimates presented in this study, wastewater data were quantitative and could accurately estimate COVID-19 prevalence. This improved correlation may also explain why wastewater has shown significantly better estimation capabilities when compared with COVID-19 hospitalizations and deaths.^2,5,49 These metrics also have lower uncertainties than reported COVID-19 cases, which are likely to be more strongly impacted by unequal distribution of testing availability, asymptomatic individuals, and individuals who do not seek testing for various reasons, including being vaccinated.

In addition to estimating community-wide prevalence with high precision, this study demonstrated that quantifying wastewater SARS-CoV-2 was a powerful method for detecting infection “hot spots” at the micro-sewershed (i.e., neighborhood) level. It should be noted that owing to the hierarchical flows between pump stations, the Northside micro-sewershed was receiving flow from the Bayfront micro-sewershed, which contains the Samaritan Pacific Communities Hospital (Figure 1). Thus, the Northside micro-sewershed wastewater results may have captured individuals who would not be linked to the Northside micro-sewershed by case reporting or random household surveillance (e.g., infected individuals at the Samaritan Pacific Communities Hospital).

Similar to the observations at the community level, wastewater SARS-CoV-2 concentrations taken from Newport micro-sewersheds correlated slightly more strongly with nasal swab positivity rates (taken during the prevalence sampling event) than with reported COVID-19 cases at the neighborhood level. Thus, wastewater SARS-CoV-2 concentrations were generally more accurate at identifying hot spots of COVID-19 infections within a community than reported clinical cases. Other studies have used similar approaches to demonstrate the utility of wastewater monitoring at the building level on college and university campuses.^56,57

WBE was also successful at identifying the dominant and minor SARS-CoV-2 variants (B.1.399/NA and B.1/NB, respectively) at both the micro-sewershed and city levels in Newport. At the city WWTP, as well as at three of the four micro-sewersheds that contained positive nasal swab samples during the random door-to-door sampling events (Big Creek, Nye Beach, and Northside), the SARS-CoV-2 variants identified in the wastewater were the same as those identified via nasal swab analyses and at similar relative proportions. The lone exception was the Bayfront micro-sewershed whose wastewater sample was dominated by the B.1.399/NA variant, whereas the nasal swab samples from that area were dominated by the B.1/NB variant. However, this micro-sewershed was home to the Samaritan Pacific Communities Hospital, which may have skewed the micro-sewershed variant information with inputs from hospitalized individuals who would not be normally associated with that micro-sewershed.

WBE data was also able to monitor the change in SARS-CoV-2 variants present in the Newport community over time. The five major SARS-CoV-2 variants identified in the WBE data from 8 June to 30 November 2020 mirrored the primary SARS-CoV-2 variants identified through clinical samples in the state of Oregon during that same time period. Interestingly, spikes observed in the WWTP influent SARS-CoV-2 concentrations often corresponded to the appearance of a new dominant variant in the wastewater sequences. This suggests that rises in viral RNA concentrations in wastewater samples may signal the appearance of a new dominant variant.

The agreement between wastewater sequence data with clinical sequence data at both the micro-sewershed and city levels supports the reliability of wastewater sequence data to accurately represent the SARS-CoV-2 variant distribution in a city or neighborhood. To our knowledge, this is the first study of variant abundance in wastewater at the neighborhood level, although others have sequenced SARS-CoV-2 in wastewater at the community level.^58–60

Although results of this study have demonstrated the power of WBE methodologies, care should be taken when comparing quantitative results across studies. A large interlaboratory study found that the quantified SARS-CoV-2 concentrations in a single wastewater sample can vary by several orders of magnitude, depending on the concentration methods (e.g., polyethylene glycol precipitation, ultrafiltration, direct extraction, or HA membrane filtration), pretreatments (e.g., removal of solids and/or pasteurization of the wastewater samples), and PCR platforms (e.g., digital or quantitative) employed.¹⁸

However, this same study also reported very low variability in replicate samples from each participating laboratory, as well as between participating labs that used the exact same methodologies, ($\text{SD}<0.2{\log }_{10}\text{\hspace{0.17em}gc}/\mathrm{L}$). This finding suggests that when the same methodologies were employed, consistent SARS-CoV-2 concentrations both within and across laboratories was achievable. Further research to optimize and standardize the wastewater surveillance process is still needed to unify the results from national and global wastewater surveillance efforts.