Year-long wastewater monitoring for SARS-CoV-2 signals in combined and separate sanitary sewers.

COVID-19 wastewater-based epidemiology has been performed in catchments of various sizes and sewer types with many short-term studies available and multi-seasonal studies emerging. The objective of this study was to compare weekly observations of SARS-CoV-2 genes in municipal wastewater across multiple seasons for different systems as a factor of sewer type (combined, separate sanitary) and system size. Sampling occurred following the first wave of SARS-CoV-2 cases in the study region (June 2020) and continued through the third wave (May 2021), the period during which clinical testing was widely available and different variants dominated clinical cases. The strongest correlations were observed between wastewater N1 concentrations and the cumulative clinical cases reported in the 2 weeks prior to wastewater sampling, followed by the week prior, new cases, and the week after wastewater sampling. Sewer type and size did not necessarily explain the strength of the correlations, indicating that other non-sewer factors may be impacting the observations. In-system sampling results for the largest system sampled are presented for 1 month. Removing wet weather days from the data sets improved even the flow-normalized correlations for the systems, potentially indicating that interpreting results during wet weather events may be more complicated than simply accounting for dilution. PRACTITIONER POINTS: SARS-CoV-2 in wastewater correlated best with total clinical cases reported in 2 weeks before wastewater sampling at the utility level. Study performed when clinical testing was widespread during the year after the first COVID-19 wave in the region. Sewer type and size did not necessarily explain correlation strength between clinical cases and wastewater-based epidemiology results. Removing wet weather days improved correlations for 3/4 utilities studied, including both separate sanitary and combined sewers.

INTRODUCTION

SARS‐CoV‐2 is shed in the feces and other fluids of COVID‐19 patients including asymptomatic individuals, making it a suitable target for wastewater‐based epidemiology (WBE). Reports of the timing and duration of fecal shedding indicated that pharyngeal swabs produced earlier positive results compared with feces (median 6.5 vs. 11 days, respectively, following symptom onset) and that shedding in feces could continue for 6–10 days after a negative pharyngeal swab test in 64% of patients in one study (Chen et al., 2020). A modeling study of reported shedding results indicated that shedding in feces could occur for 20–32 days with a median concentration in feces of 2.6 log copies/g feces (Miura et al., 2021). Many studies conducted in spring of 2020 demonstrated proof of concept for using wastewater monitoring of SARS‐CoV‐2, and examples including 2 to 4 months of data were reported for North America (Carrillo‐Reyes et al., 2021; Gerrity et al., 2021; Weidhaas et al., 2021; Wu et al., 2021), South America (Prado et al., 2021), Europe (Agrawal et al., 2021; Castiglioni et al., 2022; Medema et al., 2020; Saguti et al., 2021), Middle East (Hasan et al., 2021), and Asia (Kumar et al., 2021). Many of these studies were performed when rates of individual testing were low in many regions, worldwide one test per 1000 people in a country were reported to be available anywhere from 3/23/2020 to 12/23/2021 (with some countries still not meeting that rate) with a median date of 9/12/2020 (Ritchie et al., 2020), and thus, biomolecular detection of the virus in wastewater preceded and therefore predicted cases and hospitalizations. As testing became more widespread, wastewater monitoring has continued and grown around the world, allowing for capturing data from across seasons and waves of COVID‐19 cases (Table 1).

TABLE 1

Summary of SARS‐CoV‐2 wastewater‐based epidemiology studies that spanned multiple seasons

Authors	Location	Time period	Total number of samples	Samples/week	Duration (months)	Number of WWTP	Sewer types	Flow (MGD)	SARS‐CoV‐2 gene targets	Normalization tested
(Ai et al.,  2021 )	OH, USA	07/2020–01/2021	250	2	5	9	Both	2–130	N1, N2, E	PMMoV, crAssphage
(Feng et al.,  2021 )	WI, USA	8/30/2020–1/20/2021	418	1–5	5	12	Both	3–330	N1, N2	BCoV, PMMoV, HF183
(Li et al.,  2022 )	NV, USA	9/2020–9/2021	541	7	12	3+ in network	Separate	1.7–24	N1, N2	PMMoV
(Galani et al.,  2022 )	Greece	8/31/2020–3/21/2021		7	6	1	ND	190	N1, N2
(Koureas et al.,  2021 )	Greece	10/25/2020–04/14/2021	101		6.5	2	ND	ND	ORF1ab, S, N	Water quality
(Giraud‐Billoud et al.,  2021 )	Argentina	07/2020–1/2021	60	<1–1	6	2	ND	34–40	N1, N2
(Barrios et al.,  2021 )	Argentina	06/2020–04/2021	172		9	2+ in network	ND	ND	N1	BCoV
This study	NJ, USA	6/15/2020–5/26/2021	170	1	11	4	Both	<5–>200	N1, N2	Water quality, PMMoV

Abbreviations: BCoV, bovine corona virus; MGD, million gallons per day; ND, not described; PMMoV, pepper mild mottle virus; WWTP, wastewater treatment plant.

Studies reporting results for multi‐seasonal monitoring that occurred during periods of higher COVID‐19 testing include data from twice weekly sampling at nine wastewater treatment plants (WWTPs) in Ohio (Ai et al., 2021), once to five times weekly sampling at 12 WWTPs in Wisconsin (Feng et al., 2021), daily sampling at three WWTPs in Nevada (Li et al., 2022), and daily sampling at one WWTP in Greece (Galani et al., 2022), sampling two WWTPs in Greece (Koureas et al., 2021) and Argentina (Barrios et al., 2021; Giraud‐Billoud et al., 2021). These studies have reported correlations with reported daily cases (Ai et al., 2021; Feng et al., 2021) as well as rolling averages of 3–7 days (Ai et al., 2021), and cumulative cases for 7–20 days (Barrios et al., 2021; Koureas et al., 2021) in the sewer catchments. The correlations observed were generally consistent across WWTP; however, Feng et al. (2021) noted that one WWTP had weaker correlations with reported cases and different catchment demographics (lower income, student population) were considered a potential explanation.

While sewer systems are being used as a tool for WBE, these systems were designed to meet other goals (e.g., sanitation and hydraulics) and present multiple challenges when it comes to interpretation of results. Perhaps best studied are methods for normalizing SARS‐CoV‐2 gene copies in sewage to other system flow, fecal viruses (Ai et al., 2021; Barrios et al., 2021; Feng et al., 2021; Wu et al., 2020), fecal indicators (Feng et al., 2021), and water quality indicators (Koureas et al., 2021) to help account for dilution in these systems that collect more than feces and urine. Among the studies that spanned more than one season (i.e., at least 5 months) reported to date, a variety of correlations between wastewater virus concentrations and clinical cases have been reported. Daily to a range of 3‐ to 20‐day rolling averages of case reports were commonly tested (Ai et al., 2021; Barrios et al., 2021; Galani et al., 2022; Giraud‐Billoud et al., 2021; Koureas et al., 2021) as well as comparisons to intensive care unit admissions (Galani et al., 2022) across systems with a wide range of flow rates (Table 1).

Other factors that may impact the interpretation of sewage for system‐wide WBE could be system size/hydraulics, age, and sewer type, among others. Sewer system size was demonstrated to be sufficient to facilitate decay and therefore may impact interpretation in collection systems with longer travel times for dry weather flows with a maximum of ~1200 min (McCall et al., 2022) but was not a factor that explained the varying strength of correlations with case data across different sewer systems (Feng et al., 2021). Understanding the potential impact of sewer type is limited to select US studies that included both combined and separate sanitary sewers where no differences in correlations with cases were reported by sewer type (Ai et al., 2021; Feng et al., 2021). The potential for delay of transport due to sorption in biofilms is less clear, but accumulation of SARS‐CoV‐2 gene copies has been demonstrated in sewer biofilms (Morales Medina et al., 2022). Further advice on how to interpret data collected on wet weather days and whether that varies by system type is not clear.

The objective of this study was to compare the results of municipal wastewater monitoring across sewer types and system size. To do this, weekly wastewater samples were collected from WWTP influent sites and the SARS‐CoV‐2 N1 and N2 gene copies were quantified with qPCR. Wastewater data were compared with cases reported in the study regions. The potential impact of wet weather on sampling dates was explored as a function of sewer type and system size.

METHODS

Study site description and sewage monitoring

Sampling was performed at four WWTPs in the eastern United States, two served by separate sanitary sewers and two with sections of combined sewers. Sampling began at select locations in June or July 2020 and continued on a weekly basis through March or May 2021; details are provided in Table 2. Samples were collected at the WWTP intake weekly for all systems on Tuesdays or Wednesdays. Samples were also collected in‐system for WWTP‐D for 6 weeks (2/16/2021–3/23/2021) at three chambers along the main trunk line and from a chamber on a secondary line that fed into the main line prior to the WWTP intake. The 24‐h composite samples (1 L) were collected using automatic samplers by each partner utility. After collection, the samples were stored at 4°C and transported to the lab in a cooler with ice. Upon arrival at the lab, the wastewater samples were pasteurized in a water bath at 55°C for 1 h or in an incubator at 55°C for 90 min to inactivate viral particles (Wu et al., 2020). Pasteurization has been reported to decrease or have no impact on SARS‐CoV‐2 gene copies (as discussed by Islam et al., 2022). Immediately after pasteurization, the samples were processed for concentration and precipitation of viral particles, as described below. For quality assurance, 18 field blanks consisting of deionized water were included, one for each week during sampling as part of our dormitory COVID‐19 monitoring previously described (Fahrenfeld et al., 2021) and processed in parallel with the wastewater samples.

TABLE 2

Sampling site information including design flow (million gallons per day [MGD]), sewer type, sampling period, and estimates of industrial flows and sewer travel time

WWTP	Design flow (MGD)	Sewer type	Number of samples	Sampling period		Industrial flow	Travel time (h)
A	<10	Separate	50	June 15, 2020	May 26, 2021	negligible	2.5–3
B	<25	Separate	40	June 15, 2020	March 23, 2021	NA	NA
C	<50	Combined	44	July 16, 2020	May 26, 2021	5%	4.2–6.6
D	330	Combined	36	July 15, 2020	March 23, 2021	2.9%	4

Abbreviations: NA, not available; WWTP, wastewater treatment plant.

Concentration and precipitation of viral particles in wastewater

Concentration and precipitation of viral particles from wastewater was done following a polyethylene glycol (PEG) precipitation method (Wu et al., 2020). This method was used throughout the study period for consistency while several researchers reported method comparisons and potential improvements to this method that should be considered for future research (e.g., Kaya et al., 2022; Kevill et al., 2022). Briefly, ~200–300 ml of wastewater was filtered through a 0.22‐μm mixed cellulose esters (MCE) membrane (Millipore Sigma, St. Louis, MO, USA). The filtrate was then precipitated with PEG 8000 (8% wt/vol) and NaCl (0.3 M). The samples were shaken by hand and incubated at room temperature (20–22°C) for ~15 min until the chemicals were dissolved. The sample solution (200 ml of filtrate with PEG and NaCl) was centrifuged for 2 h at 12,000×g at 4°C immediately or after holding overnight at 4°C. Then, ~150 ml of supernatant of each falcon tube was discarded, and the remaining 50 ml was vortexed and transferred to 50 ml of falcon tube. The 50‐ml falcon tube was then centrifuged at 12,000×g for 45 min or until a pellet was visible. The supernatant was discarded and the pellet was used for total RNA extraction.

RNA extraction from wastewater pellets

Total RNA and DNA extraction of the wastewater pellets was performed using the RNeasy® PowerWater® kit (Qiagen, Germantown, MD, USA) or when supply chain issues delayed kit delivery in late 2020 (November 4 to December 9), via direct precipitation (Wu et al., 2020). When using the kit, the wastewater pellet was resuspended with 990 μl of the PM1 solution from the kit and 10 μl of 2‐mercapto‐ethanol and then transferred to the PowerWater DNA Bead tube. The rest of the protocol was completed following the manufacturer's instructions for simultaneous total RNA and DNA extraction. Briefly, for direct precipitation, wastewater pellets were resuspended in a chloroform/TRIzol reagent mixture and centrifuged to remove debris; the suspended RNA was precipitated with isopropanol, washed with ethanol, and recovered with diethyl pyrocarbonate (Wu et al., 2020). With both protocols, total RNA and DNA suspension was stored at −80°C in plastic O‐ring tubes, 100 μl of RNA free water for kit method or 50 μl for direct precipitation. Samples without detection were re‐run, and if there was still no detection, extracts were subjected to matrix spikes and re‐analyzed at a 1:1 dilution (extract:water, v:v).

Biomolecular analysis

Total RNA from the wastewater samples was analyzed by reverse transcription quantitative PCR (RT‐qPCR) to detect and quantify the SARS‐CoV‐2 N1 and N2 gene copies in the samples. RT‐qPCR was also performed for the pepper mild mottle virus (PMMoV) (Haramoto et al., 2013) to serve as an RNA extraction control because it is prevalent in sewage. The reverse transcription and PCR reactions were performed using the iTaq Universal Probes One‐step kit (Bio‐Rad, Catalog No. 1725141, Hercules, CA, USA) and the US CDC N1 and N2 primer‐probe (Table S1, IDT, Coralville, IA, USA). The PCR reaction mixture consisted of 10 μl of 2x iTaq PCR reaction mix, 0.5 μM of forward primer, 0.7 μM of reverse primer, 0.2 μM of the FAM probe, 0.25 μl of iScript reverse transcriptase, 4 μl of RNA template (sample), and RNase free ddH2O to a final volume of 20 μl. The RT‐qPCR was carried out on a Real‐Time Thermocycler (CFX96 Touch, Bio‐Rad). The one‐step PCR conditions for both N gene primer‐probe sets were the following: a reverse transcription cycle of 10 min at 50°C, a denaturing step of 10 min at 95°C, followed by 45 cycles of 10 s at 95°C and 30 s at 60°C. The one‐step PCR conditions for the PMMoV primer set included 10 min at 50°C then denaturing for 10 min at 95°C followed by 45 cycles of 10 s at 95°C and 1 min at 52°C and 1 min at 72°C. A commercial plasmid containing the complete SARS‐CoV‐2 nucleocapsid gene was used as positive control for the construction of a five‐point standard curve for both N1 and N2 (IDT, Coralville, IA). Lyophilized PMMoV (agida, Elkart, IN, USA) was extracted as described above, quantified via Nanodrop, and gene copy concentration was estimated using the equation described by Rosario et al. (2009) to serve as the PMMoV positive control for the construction of the five‐point standard curve. The range of the standard curves was 105 to 101. All standards, samples, and a no‐template control (NTC) were analyzed in triplicate (technical replicates) on each 96‐well qPCR plate. The average R 2 of the standard curves and efficiency in all reactions were 0.98 ± 0.01 and 93.3 ± 5.10% for N1, 0.98 ± 0.01 and 92.5 ± 4.8% for N2, and 0.99 ± 0.005 and 89.5 ± 4.7% for PMMoV, respectively.

The limit of quantitation (LOQ) for wastewater samples was 0.63 to 1.25 gene copies/ml as determined based on the lowest qPCR standard and the volume of wastewater filtered. The limit of detection (LOD) was the same as the LOQ, as positive control dilution lower than 101 copies/ml resulted in no detection for the qPCR. Agarose gel electrophoresis was used to confirm the presence of amplicons with the correct insert length and detection below the LOQ and LOD.

Water quality measurements

Conductivity was measured using a calibrated multimeter (Orion Star A329, Thermo Scientific, Waltham, MA, USA), and pH was measured using a calibrated Oakton pH 700 (Oakton Instruments, Vernon Hills, IL, USA). Total suspended solids (TSS) were measured in sample aliquots following the Environmental Sciences Section (ESS) Method 340.2 (Wisconsin State Lab of Hygiene, 1993). Chemical oxygen demand (COD) was analyzed according to Hach Method 8000 with Hach COD vials (20–1500 mg/L range) and a DR2700 spectrophotometer (Hach, Loveland, CO, USA).

COVID‐19 cases in the sewer catchments

COVID‐19 cases were estimated by using publicly available county‐level data (NJ COVID‐19 Dashboard) and multiplying by the percentage of the county population living in the towns served by a given WWTP. (Note that the utility catchment boundaries do not necessarily coincide with town or county borders.) US census data were used to estimate the population for the county and each town served by the WWTP. When more than one county was served by an individual WWTP, the COVID‐19 cases estimated as just described were summed for all counties served by the WWTP. Given that the study period did not include any periods of within‐state travel bans and that the towns and cities served by the WWTP studied receive wastewater from nonresidents including visitors, travelers (airports and regional rail lines pass through the study region), and workers who live in other towns and counties in this densely populated region, the residential population is not exactly representative of the population contributing to the sewage samples. Data for population movement compared with baseline during the study period for residential, workplace, and transit stations in the counties served by the WWTP studied are shown in Figure S1 (https://www.google.com/covid19/mobility/, accessed 9/28/2021). This study was reviewed by the Rutgers Institutional Review Board as secondary research and determined to be exempt following review.

Data analysis

Statistical tests were performed in R (www.r-project.org). To compare between the WWTP sampled for SARS‐CoV‐2 observations and water quality parameters, the Kruskal–Wallis tests were used followed by a post hoc pairwise t‐test with a Bonferroni correction for multiple comparisons. A two‐proportion z‐test was used to compare the proportion of samples with N1 detections to N2 detections. Spearman's correlations were tested between SARS‐CoV‐2 detection in wastewater and the number of COVID‐19 cases reported. Spearman's rhos greater than 0.8 are very strong, between 0.6 and 0.79 are considered strong, between 0.4 and 0.59 are considered moderate, and between 0.2 and 0.39 are considered weak. For the SARS‐CoV‐2 detection, tests were performed for N1 or N2 gene copies per milliliter and N1 or N2 gene copies per person per day, as calculated using Equation (1). where Q is the WWTP flowrate measured on the sampling date in milliliters per day. The catchment COVID‐19 cases were summed for the (a) 14 days prior to, (b) 7 days prior to, (c) same day, and (d) 7 days after the wastewater sampling. To understand if removing days with higher runoff improved the correlations, correlations were also tested for the 14 days prior to and the 7 days after the WW sampling, removing all data where there was any precipitation and/or snowmelt during the day of WW sampling. Climate data were downloaded from NOAA for the nearest location to the WWTP catchments (ranging 4 km north to 44.9 km southwest). Comparisons of Spearman's rhos calculated for all WWTPs were made via a Kruskal–Wallis test. Comparisons between Spearman's rhos between individual WWTPs were made via ANOVA with a post hoc Tukey honestly significant difference (HSD) test. Data normality was determined via a Shapiro test, and homogeneity of residuals via Levene's test. Correlation testing was repeated for WWTP‐D's in‐system data for the summed COVID‐19 cases in each catchment and N1 gene copies per day or milliliter for all sampling dates and dry weather dates only (i.e., dates without any reported rainfall or snowmelt during the 24‐h composite sampling). A paired Wilcoxon rank sum test was used to compare the sum of N1 gene copies per day from the WWTP‐D in‐system samples and the N1 gene copies per day observed at the WWTP intake.

RESULTS

Observations of SARS‐CoV‐2 in municipal wastewater

Spikes in concentration for the N1 and N2 genes above LOQ began in fall 2020 with the highest concentrations observed in January 2021, corresponding to the region's third wave of COVID‐19 cases (gene copies per capita per day shown in Figure 1 and gene copies per milliliter shown in Figure S2). The N1 gene peak concentrations and timing varied between the WWTPs: the highest peak concentration (N1 gene copies per milliliter) was observed as 137 for WWTP‐B on January 19, 2021, followed by 63 at WWTP‐A on January 20, 2021, 5.4 at WWTP‐C on January 20, 2021, and 3.7 at WWTP‐D on October 16, 2020. Peak N2 gene copies per milliliter of 90 were observed on January 19, 2021, at WWTP‐B, 46 on January 20, 2021, at WWTP‐A, 3.3 on March 24, 2021, at WWTP‐C, and November 18 at WWTP‐D. The highest peak concentrations of both genes were observed at WWTP‐A and WWTP‐B, the two plants receiving wastewater from separate sanitary sewers. However, there were no significant differences in paired samples across time for either N1 or N2 gene copy concentrations (copies per milliliter) nor for gene copies per capita per day between the four treatment plants studied (p = 0.34, Kruskal Wallis test). Likewise, the percentage of samples above LOQ for the N1 gene was similar between the treatment plants studied: WWTP‐A 18% (N = 9/50), WWTP‐B 17.5% (N = 7/40), WWTP‐C 22.7% (N = 10/44), and WWTP‐D 19.4% (N = 7/36). Across all treatment plants, the N1 gene was detected more frequently than the N2 gene (p = 0.009, two‐proportion z‐test; Figure S3).

FIGURE 1

(a) N1 (blue) and (b) N2 (green) gene copy concentrations normalized by wastewater treatment plant (WWTP) flow and number of people in the sewer catchment versus sampling date for each sampling location. Both y‐axes are log‐scale. Data above limit of quantitation (LOQ) are shown as circles, data below LOQ but above detection are shown as triangles, and non‐detects are shown intersecting with the respective axes. COVID‐19 cases per capita estimated for residents of the catchment are shown for the 2 weeks prior (purple crosses) to the wastewater testing. (c) N1 gene copies per capita per day versus the COVID‐19 cases per person in the catchment for the 2 weeks prior to the wastewater testing. Symbols correspond to WWTP. (d) N1 gene copies per capita per day for dry weather days versus the previous 2 weeks of cases in each catchment prior to wastewater sampling. Symbols are same as (c).

No differences were observed in PMMoV concentrations by WWTP (p = 0.57, Kruskal–Wallis). Matrix spike recoveries were 97 ± 19% (N = 18). PMMoV was quantified in all but four samples (Figure S2), which had no RNA extract left following matrix spikes such that repeating qPCR analysis with diluted extract was not possible. Notably, of the four samples without PMMoV quantification, all had the N1 gene above detection and three had N2 above detection.

During a period of extraction kit supply shortages, RNA extraction was performed by direct precipitation (N = 24 of the 170 wastewater samples in this study). Of the four samples with incomplete troubleshooting described above and therefore no quantification of PMMoV, two were from this period. A similar number of samples resulted in non‐detects with direct precipitation for N1 (33%, N = 8/24) and PMMoV (8.3%, N = 2/24) compared with samples extracted with the kit (for N1 47%, N = 69/146; for PMMoV 14%, N = 2/146, both p ≥ 0.17, prop test). For N2, a smaller proportion of samples (21%, N = 5/24) had no detection with direct precipitation compared with samples extracted with the kit (66%, N = 97/146; p = 6.3 × 10−5, prop test), potentially due to higher case rates during that period.

Correlations between WW SARS‐CoV‐2 observations and COVID‐19 cases in the sewershed

Moderate correlations were observed between N1 gene copies per capita per day and COVID‐19 cases per capita in the sewer catchments for all scenarios tested (2 weeks prior, 1 week prior, same day, and 1 week after WW sampling) across the WWTP studied (Table S2 “All,” Figures S4–S6, all rho = 0.48–0.55, all p ≤ 1.9 × 10−11). Normalizing the N1 gene copies to either flow and population or PMMoV gene copies resulted in stronger correlations than those tested for N1 gene copies per milliliter and case rates (p < 0.002, Kruskal Wallis test with post hoc pairwise t‐test). The flow and population normalization and PMMoV‐normalized correlations were similar to one another (p = 0.16, post hoc pairwise t‐test).

Testing correlations between COVID‐19 cases reported for an individual WWTP catchment (Figure 2), Spearman's rhos were similar for correlations between COVID‐19 cases reported 2 weeks prior and 1 week prior to WW sampling for all of the WWTPs sampled (p = 0.88, Tukey's HSD). The correlations between N1 observations and the COVID‐19 cases reported before WW sampling were significantly greater than those for the COVID‐19 cases reported 1 week after WW sampling across the WWTP (all p ≤ 0.02, Tukey's HSD). Therefore, wastewater was a better lagging indicator of the reported COVID‐19 cases than leading indicator across the WWTP sampled in this study (see Section 4.1 for discussion). Correlations with new cases reported the day of WW sampling were similar to the other periods tested (all p ≥ 0.09, Tukey's HSD).

FIGURE 2

Spearman's correlation coefficients (rho) for significant correlations between SARS‐CoV‐2 gene copies per volume (copies per milliliter), copies per person per day, or pepper mild mottle virus (PMMoV) normalized copies for each wastewater treatment plant (WWTP) sampled and county‐level cases for the population in the sewer catchment. COVID‐19 cases were summed for 2 weeks before the WW testing on dry days only (−2w, d), total over the 2 weeks before the WW testing (−2w), new cases (new), total over the week prior to the wastewater testing (−1w), and total over the week after the WW testing (+1w). For N1 copies per milliliter, cases were compared directly; for the N1 copies/(capita × day), COVID‐19 cases per capita in the sewage catchment were compared.

Sampling multiple WWTPs allowed for the potential to compare how well each system's WW SARS‐CoV‐2 observations corresponded with reported COVID‐19 case rates. Interestingly, the correlations observed were similar for the smallest separate system WWTP‐A and the largest combined system WWTP‐D (p = 0.16, Tukey's HSD). The weakest correlations were observed for WWTP‐C, the second largest system served by a combined sewer, or WWTP‐B, the larger of the separate sanitary sewer systems sampled (all p ≤ 0.002, Tukey's HSD).

Rainfall events were noted on some sampling dates when N1 was not observed but expected due to high COVID‐19 case rates, for example, during December 2020 and March 2021 (Figure S7). If one assumes that false negative results were observed for sampling dates with case rates of 1:10,000 (sensitivity reported based on regression of shedding models, Hewitt et al., 2022) in the 2 weeks prior to the corresponding WW sampling without N1 observations, then 18.4% of samples for WWTP‐A, 17.5% for WWTP‐B, 22.7% for WWTP‐C, and 16.7% for WWTP‐D resulted in false negative WW observations. To address this, data for dry weather N1 observations were also compared with COVID‐19 case rates in the 2 weeks prior to WW sampling, resulting in strong correlations between wastewater observations for WWTP‐D and the estimated local case rates (Figure 2, Table S2). The Spearman rhos from the dry weather correlations were significantly greater across the individual WWTP than the correlations tested that included data from wet and dry weather days (p ≤ 0.004, Tukey's HSD) (i.e., dates with any reported precipitation or snow melt during the 24‐h composite). Repeating this analysis using the cases reported in the week following WW sampling on dry weather days resulted in significant correlations for WWTP‐A but not for the other WWTP analyzed individually or together (Table S2).

In‐system sampling results

In‐system samples were collected in late winter and early spring 2021 from WWTP‐D (Figure 3). Sampling occurred in three chambers in series on the main trunk line representing 19.2%, 32.6%, and 40.2% of the flow, on average, before entering the WWTP intake (Figure 3b). Downstream of the third chamber, a second line carrying 10.8% of the average daily flow entered, combining with the third chamber outflow, before the WWTP intake. An N1 gene balance was calculated by summing the N1 gene copies per day from the third chamber on the main line and the second line and comparing this value to that observed at the WWTP intake. No significant differences were observed in paired sum of N1 gene copies per day entering the intake and observed at the intake (p = 1, Wilcoxon Rank Sum) and the relative percentage difference between the summed and observed values 86 ± 70% (ranging from 1% to 200% relative percent difference, RPD). Spearman's correlations were observed between the 2 weeks prior cases on dry weather days across the WWTP‐D samples (p = 0.0033, rho = 0.62) but not when including data from both wet and dry sampling days (p = 0.075) (Figure 3c). (Note that, given the application of Spearman's rank correlations, the same correlations would apply for the 2‐week rolling average case rates.)

FIGURE 3

(a) N1 gene copies per milliliter (black circles, left axis) for the in‐system sampling sites from the sewer for WWTP‐D. Error bars represent standard deviation of technical replicates (N = 3). Red triangles (right axis) illustrate the relative percentage difference between the sum of gene copies per day from Chambers 3 and 4 compared with the WWTP‐D influent gene copies per day. (b) Flow chart illustrating sampling chambers and the percentage of total WWTP‐D flow that pass through these chambers on sampling dates. (c) N1 gene copies per day (shapes and colors correspond to the sampling locations) versus the previous 2 weeks' sum of cases in each catchment. WWTP, wastewater treatment plant

Water quality

Water quality data were collected throughout the study period indicating some differences (Figure S8), not all of which were of practical significance. For example, pH was higher for WWTP‐C than the other plants studied (all p ≤ 0.001, pairwise t‐test); however, this difference is unlikely of practical significance as this corresponded to a difference in average pH of 0.20–0.34 pH units across the study period. Likewise, TSS was greater for WWTP‐C compared with WWTP‐A and WWTP‐B and (both p ≤ 0.035, pairwise t‐test) with no other differences observed (all p ≥ 0.14, pairwise t‐test). Conductivity was different for all WWTP pairings (all p ≤ 0.010, pairwise t‐test) except WWTP‐C and WWTP‐D (p = 1.0, pairwise t‐test). Finally, no differences were observed for COD between WWTPs (all p > 0.23, pairwise t‐test). Normalizing the N1 gene copies per milliliter or per day to the water quality parameters for correlations with the 2 weeks' or 1 week's prior cases or following week's cases resulted in similar Spearman's rhos to the non‐normalized data (Table S3).

DISCUSSION

Correlations between wastewater SARS‐CoV‐2 observations and reported COVID‐19 cases in the catchment

Wastewater SARS‐CoV‐2 concentrations observed during the second and third waves of cases in the study region, when testing was widely available but vaccine access was limited to high‐risk groups, correlated with all scenarios tested, but the strongest correlations were observed with the cases reported in the 2 weeks prior to wastewater sampling. This may be explained by observations that fecal shedding can last for 20–30 days and peak fecal concentrations for COVID‐19 patients have been reported ~10 days after symptom onset (Miura et al., 2021). This result is most similar to that reported by Feng et al. (2021) for a cross‐WWTP study in Wisconsin, USA, from late August 2020 through early January 2021, which noted either no‐lag or clinical cases peaking slightly prior to wastewater. One other US study reporting correlations for system‐level wastewater sampling campaigns noted the best correlations with daily new clinical cases in Ohio with twice weekly sampling (Ai et al., 2021) for a similar time period to the WI study. The two mid‐western studies reported strong average Spearman's rhos (OH average rho was 0.70 across 8 WWTPs, and WI average rho for per capita normalized data was 0.74 across 12 WWTPs) compared with moderate correlations observed here for New Jersey with average rho of 0.52 across four WWTPs for new cases. The other studies referenced did not specifically report testing for WW as a lagging indicator, and including such comparisons would be useful to better understand the practical significance of the differences observed in this study. This along with sampling frequency should be considered when comparing the results, as the once weekly wastewater sampling performed here may have been too infrequent to capture pre‐clinical shedding. Cross‐seasonal (i.e., >5 months) studies reporting wastewater as a leading indicator of clinical cases were performed in Greece (Galani et al., 2022) and Argentina (Giraud‐Billoud et al., 2021), and Nevada (Li et al., 2022), notably two of these studies performed daily sampling.

Limited in‐system sampling was possible in the present study but did demonstrate reasonable SARS‐CoV‐2 gene copy balance between the in‐system samples, providing proof of concept. Sub‐system sampling has been demonstrated in several other more robust studies (Baldovin et al., 2021; Barrios et al., 2021; Haak et al., 2022; Li et al., 2022) for understanding spatio‐temporal profiles.

Vaccines for COVID‐19 became available in the study region beginning mid‐December 2020 for select community members with periodic increases in eligibility among adults based on their health and job risks until all adults were deemed eligible on April 19, 2021 (Figure S9). This coincided with the third wave of COVID‐19 cases in the state, and by the end of May 2021, 50.5% of the total population in New Jersey was fully vaccinated (Figure S8). The spring 2021 lower case rates and wider availability of vaccines, along with budget limitations, all factored into the decisions of stopping the wastewater monitoring. The vaccines themselves should not result in shedding of N1 and N2 gene copies. However, the correlations reported in the present study should be considered in the context of the study period and, for example, whether the individual testing and fecal shedding would continue at similar rates as vaccination rates increased.

Sewer type, size, wet weather, and other factors

Regarding the sewer system factors (i.e., system size and sewer type) that could explain the correlations observed, correlations were similar for the smallest sampled separate sanitary and the largest sampled system with combined sewers. Therefore, as noted by others, sewer size sewer type did not necessarily explain the strength of the correlations with clinical cases observed (Ai et al., 2021; Feng et al., 2021). While one may expect the longer transport time in larger sewer systems to potentially decay SARS‐CoV‐2 gene copies, as previously demonstrated for a one‐system study (McCall et al., 2022), the travel times in the present study were shorter than those reported to result in significant decay by McCall et al. It is possible that the correlations were affected in this present study by the sampling design (once per week rather than multiple weekly sampling) and/or the factors beyond the sewer itself (e.g., demographic/geospatial variation in testing) that could impact the clinical testing data accuracy, discussed further below. The latter may be a motivation for exploring alternative clinical metrics for comparison, such as hospitalization rates or intensive care unit admissions (Galani et al., 2022), that may be less biased.

Interestingly, removing data collected on wet weather days improved the correlations for both separate and combined sewer systems in our study, beyond normalizing gene copy data for flow measured at the WWTP inlet. This improvement was not limited to the combined sewers, which may be due in part to the fact that many separate sanitary systems can have significant infiltration and inflow that could contribute to dilution (Belhadj et al., 1995), introduction of other SARS‐CoV‐2 sources (e.g., domestic animal feces given reports of positive rectal swab samples, Calvet et al., 2021), and potentially different degradation patterns. If the higher flow rates during wet weather also induce more shear stress on sewer biofilms and sediment, they could result in re‐introduction of previously shed virus found in these matrices (Fahrenfeld et al., 2021; Morales Medina et al., 2022). But more studies would be needed to demonstrate the potential impact of such a phenomenon beyond simply demonstration of SARS‐CoV‐2 gene copies in these matrices. Two other studies that sampled both combined and separate sanitary sewers did not report differences in correlations by sewer type (Ai et al., 2021; Feng et al., 2021), although one did note a general lack of major precipitation events on sampling dates during the study period (Ai et al., 2021). To account for the potential dilution during wet weather events, others have also normalized to flow and/or other fecal makers to account for dilution. The highest influent flow rates were not necessarily observed on the wet weather days, which may have been masked by industrial flows, lack of large rainfall events near sampling events (e.g., rainfall >2 in, > 50.8 mm), or seasonal variation in the influent flows in the catchments. Normalizing to fecal marker did not always improve the correlations in the present study, as described below.

Normalization to other marker genes was performed for PMMoV in the present study. PMMoV was not quantified in four samples, resulting in a minor loss of data for testing correlations. Using all the WWTP data and normalizing to PMMoV resulted in similar correlations to the per capita per day normalized N1 gene copy data, making the utility of PMMoV normalization less clear. A concurrent study from our lab on the building level found that normalization to the PMMoV did not improve correlations (Fahrenfeld et al., 2021), and others have reported that correlations did not improve when correcting for a spiked control or different fecal markers (i.e., PMMoV and HF183, Ai et al., 2021; Feng et al., 2021). Normalizing to water quality parameters (i.e., pH, TSS, and COD) did not improve the correlations observed here. In contrast, normalizing SARS‐CoV‐2 concentrations for COD for data from two WWTPs, but not biochemical oxygen demand (BOD) or TSS, was reported to improve correlations in a case study from Greece (Koureas et al., 2021).

Beyond sewer systems factors and testing availability, there may be multiple social and virus factors that also impact the interpretation of wastewater‐based epidemiology results. The correlations tested here against reported cases may be impacted if there is sufficient testing hesitancy among the target population, as has been reported to be associated with economics, race, and immigration status, all of which may vary geospatially (Egelko et al., 2020). One study of nine sewer systems noted that one system that performed worse had significantly different population demographics consisting of primarily university students compared with the other systems studied (Ai et al., 2021). Whether the shift in dominant variants affected the results presented here is difficult to access with the presently available data, as the state of New Jersey began reporting variant estimates mid‐way through the study period (January 2021).

CONCLUSIONS

The present work provides a cross‐seasonal study of SARS‐CoV‐2 wastewater‐based epidemiology during which testing was widely available and different dominant variants were reported in the clinical population. Under these conditions, weekly wastewater concentrations correlated best with clinical data from the 2 weeks and 1 week prior to wastewater collection or new cases rather than cases reported 1 week after. The strength of these correlations was similar between the sewers of different sizes and types. Future studies should take advantage of the method development performed by other teams during the study period that can result in higher virus yield (e.g., Kaya et al., 2022; Kevill et al., 2022) as well as in comparison with studies reporting higher sampling frequency. Nonetheless, removing wet weather dates improved the correlations with clinical cases for three of the full systems studied and the short‐term in‐system study, even when already normalizing data for flow, which may indicate that interpretation of weekly wastewater epidemiology data during wet weather events may be more complicated than simply accounting for dilution.

AUTHOR CONTRIBUTIONS

Nicole L. Fahrenfeld: Conceptualization; data curation; formal analysis; funding acquisition; investigation; methodology; project administration; supervision; validation; visualization. William Morales Medina: Data curation; formal analysis; funding acquisition; investigation; methodology. Stephanie D'Elia: Data curation; formal analysis; investigation; methodology; validation. Aishwarya Deshpande: Investigation; methodology. Genevieve Ehasz: Investigation; methodology.

Table S1 qPCR primers and calibration curve data

Table S2 Spearman correlation coefficients (rho) and p‐values for SARS‐CoV‐2 gene copies per volume (copies per mL), copies per capita per day for each wastewater treatment plant (WWTP) sampled, or PMMoV gene copy normalized N1 gene copies and county level cases for the population in the sewer catchment. COVID‐19 cases were summed for two weeks before the WW testing on dry days only (−2w, d), two weeks before the WW testing (−2w), the week before the wastewater testing (−1w), the day of WW testing (new), the week after the WW testing (+1w), and the week after WW testing on dry days only (+1w, d). For N1 copies/mL cases were compared directly, for the N1 copies/capita*day, COVID‐19 cases per person in the sewage catchment were compared.

Table S3 Spearman correlation coefficients (rho) and p‐values for SARS‐CoV‐2 gene copies per volume (copies per mL) or copies per capita per day normalized to water quality data across the WWTP sampled and county level cases for the population in the sewer catchment. COVID‐19 cases were summed for two weeks before the WW testing (−2w), the week before the wastewater testing (−1w), or the week after the WW testing (+1w). For N1 copies/mL cases were compared directly, for the N1 copies/capita*day, COVID‐19 cases per person in the sewage catchment were compared.

Figure S1 Changes from baseline in movement for the counties (represented by different colored lines) served by the WWTP (A‐D) studied here for (a) residential, (b) workplaces, and (c) transit stations. Some counties were served by more than one WWTP, as indicated by there being two letters in the graph title frame. The vertical red dotted line indicates when the first wastewater sampling for this study began. Data from https://www.google.com/covid19/mobility/ accessed 9/28/2021.

Figure S2 Log10 (a) N1 gene copies/mL, (b) N2 gene copies/mL, and PMMoV gene copies/mL versus time for each WWTP sampled. Dotted red line is for the maximum LOQ (0.63–1.25 gene copies/mL). Error bars represent standard deviation of technical replicates (N = 3).

Figure S3 Log10 transformed concentration of N2 gene copies versus N1 gene copies with colors corresponding to different sampling sites. Red lines represent the maximum LOQ. Non‐detects are shown intersecting with the respective axes.

Figure S4 COVID‐19 reported cases per capita in the sewer catchment across the study period. Purple crosses represent cumulative cases two weeks prior to wastewater sampling, red triangles represent cumulative cases one week prior to wastewater sampling, and orange circles represent cases in the week after wastewater sampling.

Figure S5 Flow and population normalized concentrations of N1 gene copies versus the positive COVID‐19 cases per capita estimated for the (a) dry weather days two weeks before, (b) two weeks before (c) the week before, (d) new cases, and (e) the week after wastewater collection. Note, both axes are in log10 units.

Figure S6 PMMoV gene copy normalized concentrations of N1 gene copies versus the positive COVID‐19 cases per capita estimated for the (a) dry weather days two weeks before, (b) two weeks before (c) the week before, (d) new cases, and (e) the week after wastewater collection. Note, both axes are in log10 units.

Figure S7 N1 gene copy concentration (black circles, log‐scale, left axis) for sampling dates. Red dotted line represents the quantitation limit. Rainfall (blue triangles) for the nearest rain gauge on the sampling date.

Figure S8 Water quantity and quality data across time for the four wastewater treatment plants (a) percentage of design flow observed on the sampling date at the WWTP intake, (b) pH, (c) total suspended solids (TSS), (d) chemical oxygen demand (COD), and (e) conductivity.

Figure S9 The number of people (left y‐axis) receiving their first dose of vaccine for COVID‐19 (orange bars) and completing vaccination (blue bars) versus time for the study period in the state of NJ. The cumulative percent the total persons in NJ who were fully vaccinated (black line, right y‐axis) across time is overlaid. Data from US CDC and NJ COVID‐19 Dashboard.

Click here for additional data file.