SARS-CoV-2 in wastewater settled solids is associated with COVID-19 cases in a large urban sewershed

Sample collection and storage

Influent and primary settled solids were collected from two publicly owned treatment works (POTWs), hereafter referred to as POTW A (Palo Alto Regional Water Quality Control Plant) and POTW B (San José-Santa Clara Regional Wastewater Facility), both located in Santa Clara County, California, USA. POTW A and B serve populations of 0.2 million and 1.5 million, respectively with permitted flows of 39 and 167 million gallons per day. POTW B adds FeCl₃ to its waste stream (10 mg/l) prior to the headworks for odor control. The residence time of sludge in the primary clarifier is estimated to be between 2-14 hours for POTW A and 1-2 hours for POTW B.

From both plants, 50 ml of influent was collected at the headworks as 24-hour flow-weighted composite samples and 50 ml of solids samples were collected from the primary settling tank. Grab solids samples were collected from POTW A and 24-hour composite samples from POTW B. Samples were collected in 10% HCl acid-washed plastic containers. Storage conditions are described in the Supporting Information (SI).

Paired influent and solids samples were collected at both POTWs on the following dates: 22 March, 25 March, 29 March, 1 April, and 15 April 2020. Additional dates were included for POTW B (23 March and 30 March) (Table 1). Dates were chosen to span a high prevalence period in the initial phase of the pandemic. These samples are hereafter referred to as “method evaluation samples”.

A longitudinal collection of primary settled solids was obtained from POTW B. Eighty-nine samples were collected daily from March 16 to May 31, and three times a week from June 2 to July 12, using the same collection techniques as the method evaluation samples. Hereafter, these are referred to as “longitudinal samples”.

Influent pre-analytical processing

Immediately prior to analysis, influent samples were transferred to 4°C until thawed (range: 12-48 hours). Between 43 and 45 mL of influent was centrifuged at 24,000xg for 15 minutes at 4°C to pellet solids. 40 to 42 mL of the resultant clarified supernatant was decanted and viruses were concentrated from the supernatant using a PEG precipitation method²¹. This method was selected after preliminary head-to-head testing of PEG precipitation and centrifugal ultrafiltration to concentrate viruses from influent yielded comparable recovery results for N1 and N2 (Figure S1), and because we anticipated there would be fewer pandemic-related supply chain issues with the PEG method. MHV strain A59 and PMMoV were used as an exogenous and endogenous recovery controls, respectively. PMMoV also served as a fecal-strength control. An attenuated vaccine strain of bovine coronavirus (BCoV) was used as a nucleic-acid extraction positive control; it was spiked into the viral concentrate before it was subjected to nucleic-acid extraction using the Qiagen AllPrep PowerViral DNA/RNA kit and further purification using Zymo OneStep PCR Inhibitor Removal Columns (Zymo Research, Irvine, CA). Further details are in the SI.

Solids pre-analytical processing

Immediately prior to analysis, solid samples were transferred to 4°C until thawed (range: 12-48 hours). Thawed solids samples (35-53 ml) were centrifuged at 24,000xg at 4°C for 30 minutes, and the supernatant was decanted. Percent solids was measured for concentrated solid samples by placing sample aliquots in pre-weighed aluminum weigh dishes and weighing the sample before and after drying at 105°C for 24 hours.

The method evaluation solids samples were divided to test two different extraction processes: for the first, 0.2-0.4 g of sample was aliquoted for extraction with the Qiagen AllPrep PowerFecal DNA/RNA kit (“powerfecal” kit), and for the second, 1.8-2.8 g of sample was aliquoted for extraction with the RNeasy PowerSoil Total RNA Kit (“RNeasy” kit). Aliquots were spiked with BCoV as an extraction recovery control. Two biological replicates were completed per sample (2 powerfecal and 2 RNeasy extractions per POTW per time point). Nucleic-acids were further purified using the Zymo OneStep PCR Inhibitor Removal Columns (Zymo Research, Irvine, CA). Further details are in the SI.

For the longitudinal solids samples, RNA was extracted using the RNeasy kit and purified with the Zymo inhibitor removal kit using the same protocol described above with two modifications ²⁴; (1) after step 9 in the kit’s instructions, samples were stored overnight at -20°C, and (2) positive pressure was used per manufacturer’s instructions to aid flow through columns. 7 of the 89 samples were extracted in duplicate and extraction blanks were run with each sample batch. Samples were randomized and blinded to laboratory technicians prior to analysis.

RT-QPCR

For method evaluation samples, PMMoV, MHV, BCoV, SARS-CoV-2 N1 and N2 were quantified using two-step RT-QPCR. Preliminary experiments suggested that the RT step was inhibited, so dilutions of the extract were used as template for the RT step (see SI for details). For influent, undiluted and 1:10 diluted RNA was used as template; for solids, undiluted, 1:10, and 1:50 diluted RNA were used as template. Each template was run in triplicate, and triplicate reactions were pooled prior to QPCR.

Resultant cDNA was used as template in QPCR assays targeting PMMoV, MHV, BCoV, SARS-CoV-2 N1 and N2 (primers and probes in Table S1) using an ABI Step-One Plus Instrument (Thermo Fisher Scientific, Waltham, MA). All templates were run in duplicate. Standard curves were run in duplicate on each plate using cDNA at concentrations from 3 to 3 × 10⁵ copies / reaction; duplicate NTCs were included on each plate. Detailed methods for generation of cDNA standards are in the SI. Concentrations per reaction were converted to copies per volume of influent or per g dry solids weight using dimensional analysis. Samples were considered “below the limit of quantification” (BLOQ) if the Cq value of the sample was higher than the Cq value of the lowest cDNA standard, but less than 40, or if the targets in both technical replicates did not amplify within 40 cycles. Concentrations of unknown samples were calculated using plate-specific standard curves. Additional details are in the SI²⁵.

ddPCR

RNA from the method evaluation samples was used as template in one-step digital droplet (dd)RT-PCR for N1 and N2 SARS-CoV-2 targets using BioRad SARS-CoV-2 Droplet Digital PCR kits (cat.no.: 12013743; Hercules, CA) and a BioRad QX200 AutoDG Droplet Digital PCR system. Technical duplicates were run (2 wells); wells were merged for data analysis. In order to test for inhibition, each template was run at two dilutions: undiluted and 1:10 dilution. A positive control consisting of SARS-CoV-2 RNA isolated from a nasopharynx swab of a high-titer patient from Stanford Hospital and NTCs were run on each plate in duplicate. QuantaSoft and QuantaSoft Analysis Pro (BioRad) were used to manually threshold and export the data. The required number of droplets for a sample with merged duplicate wells was at least 10,000. In order for a sample to be scored as positive, 3 or more positive droplets were required in the merged wells. Concentrations per reaction were converted to copies per volume of influent or per g dry weight using dimensional analysis.

A two-step ddRT-PCR assay for N1 was also trialed using a subset of the method evaluation samples. RNA was stored at -80°c for up to 5 days between running the one-step and two-step assays, to minimize the potential for RNA degradation during storage (see SI).

For the longitudinal samples, one-step ddRT-PCR duplex assay for N1 and N2 and a duplex assay for PMMoV and BCoV²⁶ were used. N1 and N2 were quantified in undiluted and 1:10 diluted template; BCoV and PMMoV were quantified in 1:10 and 1:1000 diluted template. The N1/N2 duplex assay was conducted in triplicate wells and the BCoV/PMMoV duplex assay was conducted in duplicate wells. The replicate wells were merged for data analysis. Positive N1/N2 and negative controls, and threshold setting for the longitudinal samples were conducted as described for the method evaluation samples. Positive controls for BCoV and PMMoV were a direct extraction of reconstituted BCoV vaccine diluted to 10⁶ cp/mL and for PMMoV a synthetic DNA ultramer. In addition, a pooled matrix control was run on every N1/N2 duplex assay plate. This was created by adding equal-volume aliquots of each sample on the plate to one tube and spiking this pooled matrix with clinical SARS-Cov-2 RNA. Additional details of dd-RTPCR²⁷ can be found in the SI and Table S7.

COVID-19 epidemiology data

The number of new clinically confirmed COVID-19 cases within the sewershed of POTW B were obtained using geo-referenced case data and POTW-provided shape file of their service area in ARC-GIS. The date stamp is that of specimen collection. As data were handled and provided by the county, no IRB approval was needed.

Statistics

Statistics were computed using Microsoft Excel and Rstudio (version 1.1.463). Paired and unpaired t-tests were used to compare groups after confirming data were normally distributed. For calculating recoveries, and making numerical comparisons across samples using the high copy number targets measured by RT-QPCR (PMMoV, MHV, and BCoV), the results from reactions that yielded the highest measured concentration were used to balance loss of signal against effects of inhibitory substances. For ddRT-PCR data, “total error” from merged wells is reported in the form of standard deviations. Pearson correlation coefficients assessed association between different measurements among samples.

We used a resampling-based strategy to evaluate whether the observed measurements N1 and N2 with and without PMMoV normalization were predictive of clinically confirmed cases. We applied first differencing to account for autocorrelation, resulting in the following regression specification (using N1 as the example): Cases_t − Cases_t-k ~ N1_t – N1_t-k where the subscript t indicates the day of measurement and k is the number of days used for differencing, in this case 7; k = 14 was also used as a sensitivity test. If a wastewater data point was missing at day t-k, then the value for the closest date was imputed. In order to account for the technical variability of the measurements, the observed top and bottom confidence bounds were used to resample wastewater concentrations uniformly. For samples below the detection limit, we used uniform sampling between 0 and the lower detection limit (estimated at 40 copies/g dry weight). In addition to raw case counts, 7-d smoothed case data were also used to account for the day-to-day variation in case data. A downsampling analyses was conducted to investigate the frequency of sample collection needed to observe associations between case counts and wastewater data. We downsampling to fortnightly, weekly, and twice-per-week frequencies using a random observation within each appropriate period and examined as described above except k for biweekly was are 2 or 4 half-weeks, and for weekly and fortnightly the k was 1 or 2 (weeks/fortnights). Empirical p-values were determined with 1,000 resamplings were reported.

Analytical controls and assay performance

NTCs and extraction blanks were negative for all targets in all analytical methods. QPCR assay efficiencies ranged from 78-105% and standard curve R² ranged from 0.98 to 0.99 (Table S2; Figure S2).

Method Evaluation: Comparison of RT-QPCR vs ddPCR assays on influent and settled solids

With RT-QPCR, N1 and N2 were not detected in any influent sample, even when RNA template was diluted 1:10 and cDNA was diluted 1:10, except for a single influent from POTW B (Table S3) that resulted in detection below the limit of quantification. When one-step ddRT-PCR was applied to the same undiluted and 1:10 diluted RNA templates, one sample tested positive for N1 in POTW A, three tested positive for N1 at POTW B, and one tested positive for N2 in POTW B (Table S4). In all but one of the samples that resulted in positive measurements, the undiluted template yielded a non-detect, suggesting inhibition of the reactions in the undiluted template.

For solids samples extracted with powerfecal kits and analyzed with RT-QPCR, SARS-CoV-2 targets were detected in at least one biological replicate from 1 of 5 samples from POTW A and 4 of 7 samples from POTW B (Table S3). Detection was inconsistent across biological and technical replicates, and there was no case where both biological or technical replicates were positive. When analyzed with ddRT-PCR, the SARS-CoV-2 targets were detected and quantified in 0 out of 5 POTW-A solids samples and 7 out of 7 POTW B samples (Table S4). ddRT-PCR results were more consistent across biological replicates than RT-QPCR (Figure 1; Table S4).

• Download figure
• Open in new tab

Figure 1.

N1 and N2 measured at POTW B using ddRT-PCR. Top row: Concentrations in influent in units of copies (cp) per ml. Middle row: Concentrations in solids extracted using the powerfecal kit in units of cp per g dry weight. Bottom row: Concentrations in solids extracted using the RNeasy kit in units of cp per g dry weight. Open symbols indicate those where no target was detected (less than 3 positive droplets). Error bars represent standard deviations as represented by the “total error” which includes Poisson error as well as differences among merged wells. Results for undiluted and 1:10 diluted template are provided for the solids; results for influent are for 1:10 diluted template for N1 and undiluted template for N2 – results for dilutions not shown are all non-detects.

The improved performance of ddRT-PCR compared to RT-QPCR for both influent and solids samples was unlikely due to differences in template amounts. The virtual influent volumes in RT-QPCR and ddRT-PCR assays were 2 mL and 8.4 mL, respectively (when template not diluted); and for solids, the virtual wet mass of solids assays in the RT-QPCR and ddRT-PCR reactions were 20 mg and 50 mg, respectively (when template not diluted); thus theoretical lower detection limits are similar across analytical approaches (Table S5). There are several possible explanations for the low and inconsistent detection of SARS-CoV-2 targets in influent and solids with RT-QPCR, including assay inhibition, and the absence or low occurrence of viral RNA in the samples. We studied RT-QPCR inhibition in more depth using endogenous PMMoV and spiked MHV and BCoV (Figure S3). The results demonstrate inhibition for the RT-QPCR assays of endogenous and exogenous recovery targets, despite efforts to alleviate it using the Zymo OneStep Inhibitor Removal columns and dilution of template. For influent, PMMoV and MHV concentrations were generally highest when RNA template was diluted 1:10 prior to RT, compared to results from undiluted and 1:10 diluted cDNA generated from undiluted RNA template (Figure S3). This suggests the RT step, and not the PCR step, was most inhibited for the PMMoV and MHV RT-QPCR assays. For solids, extract dilutions of 1:50 appeared to alleviate inhibition for RNA obtained from POTW A (Figure S3; Table S6). For solids from POTW B, however, inhibition may not have not alleviated even after a 1:50 dilution.

Based on these results comparing the performance of RT-qPCR and ddRT-PCR, and the improved sensitivity of the ddRT-PCR, we selected the ddRT-PCR methods for further SARS-CoV-2 analysis.

Method Evaluation: Comparison of RNeasy and powerfecal extraction kits and 1-step vs. 2-step ddPCR for settled solids

We used the SARS-CoV-2 ddRT-PCR assays and the POTW B solids to compare powerfecal and Rneasy solids extraction kits. The virtual wet mass of solids assayed in the duplicate merged wells was approximately 50 mg (powerfecal) and 250 mg (RNeasy) when the templates were not diluted (theoretical detection limits in Table S5). Samples treated with the powerfecal kit were more commonly positive for SARS-CoV-2 targets in undiluted template than in 1:10 diluted template (Figure 1). By contrast, samples treated with the RNeasy kit were positive for SARS-CoV-2 at both dilutions, but detection was more consistent with the 1:10 diluted templates (Figure 1). The SARS-CoV-2 target concentrations were similar across the biological replicates and the null hypothesis that paired replicates have the same concentration was not rejected (p>0.05 for all 4 comparisons (2 targets × 2 RNA extraction kits), paired t-test). The concentrations of N1 and N2 measured in the same sample / biological replicate were not different between the two kits (paired t-test, p>0.05). Although the powerfecal kit samples less material than the RNeasy kit, these results suggest that some of the RNeasy kit extracts require a 1:10 dilution, thus negating the advantages of the extra sample mass.

Whereas the N1 target was detected in every solids sample from POTW B when the 1-step ddRT-PCR assay was employed, it was not detected in any of the RNeasy extracts when a 2-step ddRT-PCR assay was conducted.

Method Evaluation: Comparison of SARS-CoV-2 RNA concentrations measured in influent and primary sludge

We compared SARS-CoV-2 measurements made in influent and solids samples collected on the same days at POTW B (Figure 1). Solids results obtained using the powerfecal kits are used for the comparisons. At POTW B, the N1 target was detected in 3 out of 7 influent samples and 6 out of 7 solids samples and the N2 target was detected in 1 out of 7 influent samples and 7 out of 7 solids samples. The N1 target was quantified in both influent and primary solids on 3 dates (3/29, 4/1, 4/15) and the N2 target was quantified in both influent and primary solids on 1 date (3/30). The ratio of N1 solids concentrations (cp/g; mean of biological replicates) to influent concentrations (cp/mL) on the 3 dates was 720, 620, and 350 ml/g respectively. For the N2 target, the ratio was 3100 ml/g. These results suggest that on a per mass basis, these targets are present at between ~100 and ~1000 higher concentrations in solids.

Method Evaluation: Recovery of viral targets through preanalytical processing

Recovery measurements help characterize preanalytical method consistency across samples. For the influent method, we measured the recoveries of endogenous PMMoV and spiked MHV through the PEG concentration method (Figure S4). Median PMMoV recovery was 21% (range: 8% to 30%) and median MHV recovery was 7% (range: 2% to 16%). MHV recovery was lower than that for PMMoV by, on average, 13% (paired t-test, p<0.05). We used BoCoV to assess recovery specifically through nucleic-acid extraction and purification, and through inhibitor removal steps. Recovery of BCoV through these steps ranged from 3% to 48% (median= 26%). The PEG concentration step and nucleic acid extraction/purifications step are conducted in series. Assuming that the BCoV extraction recoveries apply to all viral targets, our results suggest that viral recoveries are between 0.1-7% (median 1%, based on MHV) or 1-8% (median 7%, based on PMMoV) of actual influent concentrations.

There were no significant differences in recoveries of PMMoV, MHV, and BCoV in influent between the two POTWs (t-test, p>0.05), or in samples that were positive or negative for N1 or N2 (t-test, p>0.05). In fact, the influent sample with one of the lowest MHV and BCoV recoveries (POTW B on 3/29/20) was one of the few influent samples positive for N1.

PMMoV may serve not only as a process control, but also a means of assessing the “fecal strength” of sewage. PMMoV measured in the viral concentrate varied between 7.7 × 10² to 7.0 × 10³ copies/ml of influent across plants (Figure S4). There was no significant difference in PMMoV influent concentrations between plants (p>0.05).

Method Evaluation: Endogenous and spiked virus recoveries in solids from method evaluation study

There was no viral concentration steps in the solid methods and BCoV served as an extraction recovery control. BCoV recovery varied from 0% to 33% (median = 16%) using the powerfecal kit and 0% to 48% (median = 22%) using the RNeasy kit (Figures S5 and S6). The 0% recovery values were obtained in POTW B solids collected on 15 April and treated with the RNeasy kit and POTW B solids collected on 25 March and treated with the powerfecal kit. Interestingly, PMMoV and SARS-CoV-2 targets were measured in the same extracts with 0% BCoV recovery at concentrations not largely different from other dates. The measured BCoV recoveries were not different across plants or across extraction kits (t-test, p>0.05). It is possible, however, that the BCoV recoveries in the POTW B solids are underestimates because we are not certain that the RT-QPCR assay conducted at the highest extract dilution was free of inhibition.

PMMoV in solids from POTW A and POTW B varied between 5×10⁴ and 7×10⁷ copies/g dry weight (median 2 × 10⁶ copies/g) when extracted with powerfecal kit (Figure S6). Only samples from POTW B were measured with the RNeasy extraction kit and the resulting concentrations ranged from 5×10⁵ to 3×10⁶ copies/g dry weight (median 1 × 10⁶ copies/g; Figure S5). These concentrations are twice as high, on average, than those measured using the powerfecal templates (paired t-test, p<0.05). The PMMoV concentrations may be underestimates due to unresolved assay inhibition. Concentrations of PMMoV were not different across biological replicates regardless of which extraction kit was used (paired t-test, p>0.05). PMMoV concentrations measured in solids at POTW A were not different from those measured at POTW B (t-test, p>0.05, data from powerfecal extracts).

N1 and N2 concentrations in solids are not correlated to BCoV recovery or PMMoV concentrations regardless of kit used (p>0.05) (Figures S5 and S6).

Longitudinal data from POTW B

The method performance experiments demonstrated that the primary solids method resulted in more positive and consistent SARS-CoV-2 signals than the influent method. Furthermore, we found that the 1-step ddRTPCR assay performed on the solids extracts resulted in more positive results than the 2-step RT-qPCR and 2-step ddRTPCR assays. We therefore focused our longitudinal analysis on primary settled solids samples from POTW B using the RNeasy kit and ddRT-PCR for measuring all targets (Figure 2).

• Download figure
• Open in new tab

Figure 2.

Top panel: Confirmed new cases of COVID-19 in the sewershed of POTW B (black) and 7-d smoothed new cases (red). Middle and bottom panels: Concentrations of N1 and N2 measured in solids (copy per g dry weight). Error bars are standard deviations as total errors from the ddPCR machine. Open symbols are non-detects plotted at 0. The theoretical lowest concentration measurable is ~40 copies/g dry weight.

Concentrations of N1 and N2 in each sample were derived from the reactions that used undiluted RNA template unless those that used 1:10 diluted template produced a concentration higher than the upper standard deviation for the undiluted template (10/96 and 6/96 samples required 1:10 dilution for N1 and N2, respectively).

N1 and N2 were detected in 79 and 68 of the 96 samples. When detected, average N1 and N2 values were 870 cp/g and 730 cp/g with maximums of 3330 cp/g and 4250 cp/g, respectively. Replicate extractions generally yielded concentrations that were similar.

Median PMMoV was 2×10⁷ cp/g (interquartile range 1 ×10⁷ cp/g to 4 ×10⁷ cp/g) (Figure S7) suggesting fairly constant fecal strength of the solids. We normalized measured concentrations of N1 and N2 to obtain fecal strength normalized target concentrations to include in the statistical analyses. The recovery of BcoV was lower in the longitudinal samples than it was in the method evaluation samples, with a median recovery was 3% (interquartile range 1 to 5%) (Figure SA). We omitted the N1 and N2 data from our dataset for the one sample where 0% BCoV recovery was measured (16 March). Recovery was higher than 10% on 8 days; N1 and N2 concentrations were not significantly different on those days compared to the others (p>0.05).

The concentrations of N1 and N2 tracked the reporting of new cases in the sewershed (based on specimen date) when accounting for autocorrelation and technical errors associated with the wastewater measurements (using raw case counts, 0.024 cases / copies /g for both N1 and N2, p = 0 and 0.005, respectively). Results were not sensitive to choice of k, smoothing of case data, or normalization of N1 and N2 by PMMoV to account for changing sample-to-sample differences in fecal strength or recovery (Table S8). Downsampling wastewater data collection and analysis frequency to down sampling to twice-per-week yields significant associations between case counts and wastewater concentrations; however, downsampled to fortnightly or once per week, associations are no longer significant (Table S8).