Wastewater monitoring of SARS-CoV-2 RNA at K-12 schools: comparison to pooled clinical testing data

Sample collection

Wastewater samples were collected from two schools in California where the superintendent of the school district granted permission to collect samples and informed their community members of the study. School A is a middle school with approximately 1,000 individuals including 900 students. School B is a lower elementary school for students in kindergarten to 2nd grade with approximately 400 individuals including 320 students. Both schools are located in the same county, within 5 km of each other. The two schools were chosen based on the availability of access points to sewage from the majority of the school. Wastewater was collected through access points that allowed us to sample before wastewater from the school buildings mixed with wastewater from other parts of the community (the access points were referred to as “sewage cleanout”); we estimate that wastewater spent < 30 min in the sewage pipes between where it entered the system and the sampling location. School A wastewater samples included wastewater from all buildings except for the gym and the music room, and School B wastewater samples included wastewater from the entire school.

Composite wastewater samples were collected using an autosampler (Teledyne ISCO, NE) programmed to collect 100 mL of wastewater every 5 min between 7:00 and 15:00 three times per week (Tuesdays, Wednesdays, and Fridays). This time period was selected to include students’ drop-off and pick-up times, which were around 8:00 and 14:00 respectively for both schools. Sampling was conducted for eight weeks from April 2022 to June 2022. Samples were retrieved by a field technician at the end of each day after samples were composited into a daily sample by the autosampler. The total volume for each sample varied according to the performance of the autosampler and ranged from 0 to 9.6 L. At the end of every sampling event, technicians decanted ∼8.5 L from the unmixed sample. The remaining 1 L of sample with a high solids content was immediately transported to the laboratory on ice. Once at the laboratory, the sample was stored at 4 °C and analyzed within a week. Previous work indicated minimal degradation of viral RNA RT-PCR targets during this time frame (Roldan-Hernandez et al., 2022; Yang et al., 2022). The autosamplers malfunctioned twice and no sample was available at School A on April 19 and at School B on April 22. During the last week of this study, there was a national holiday on Monday; therefore samples were taken on Wednesday, Thursday, and Friday during this week.

Pre-analytical processing

Wastewater samples were separated into a solid fraction and a liquid fraction by settling the sample in an Imhoff cone for 1 h at 4 °C (United States Environmental Protection Agency, 1999). The solid fraction was removed from the cone by decanting the liquid fraction into a sterile container and transferring the solid fraction to a different sterile container. Samples were processed using a modified version of Wolfe et al (2021b). The solid samples were dewatered by centrifuging at 24,000 ×g for 30 min at 4 °C and decanting the resulting supernatant. Solids included visible debris such as toilet paper. 0.225 g of solids were resuspended in DNA/RNA shield (Zymo Research, Irvine, CA, USA) so that the final concentration was 75 mg/mL, a concentration previously shown to minimize inhibition (Huisman et al., 2022). The DNA/RNA shield was spiked with bovine coronavirus (BCoV, Calf-guard Cattle Vaccine, PBS Animal Health, OH) at a concentration of 500,000 copies/mL. BCoV was used as a spiked-in internal control since the wastewater is generated solely by the school population without an agricultural source and we do not expect BCoV to be naturally present in the wastewater. The resuspended sample was stored at 4 °C overnight until nucleic acid extraction. If there were any leftover solids, dry weight of the dewatered solids was determined by drying the sample at 105 °C for 24 h.

Viruses from 45 mL of the previously separated liquid fraction were concentrated using electronegative filtration adapted from previous studies (LaTurner et al., 2021; Graham, Anderson & Boehm, 2021). The samples were centrifuged at 4,100 ×g for 10 min at 4 °C to eliminate any larger solids to prevent clogging. Sterile membranes (47 mm diameter; 0.45 µm pore size; Millipore HAWG047S6) were placed on sterile plastic filter funnels, wet with 1 mL of 1 × phosphate buffered saline solution (PBS buffer, Gibco, NY, USA) and the supernatant was transferred to the filter funnels. The sample was mixed with 1 mL of 1.25M MgCl₂ and allowed to sit for 5 min. The sample was then vacuum filtered (approximate vacuum pressure of 17 kPa) through the sterile membrane. The vacuum was turned off to relieve the pressure and then 0.2 mL of DNA/RNA shield was placed onto the membrane, allowed to sit for 5 min, then aspirated from the filter using vacuum. The membrane was folded and placed in a bead beating tube included in the Qiagen AllPrep PowerViral DNA/RNA Kit (Qiagen, Hilden, Germany) for extraction. 0.2 mL of DNA/RNA shield spiked with BCoV to a concentration of 500,000 copies/mL was added as a spiked-in internal control. Samples were stored at 4 °C overnight until nucleic acid extraction.

Nucleic acid extraction

For solid samples, 5/32” stainless steel grinder balls (OPS Diagnostics, NJ) were added to the suspension containing the solids and shaken with a bead beater (FastPrep-24 Tissue and Cell Homogenizer, MP Biomedicals, CA, USA) to homogenize the sample. The homogenized solution was centrifuged at 5,250 ×g for 5 min. Nucleic acids were extracted from the supernatant using the PowerViral Kit; manufacturer’s instructions were followed starting after the bead beating step. For liquid samples, membrane filters already added to the bead beating tubes in PowerViral Kit were processed as directed by the manufacturer. Qiacube (Qiagen, Hilden, Germany) was used to extract nucleic acids for both solid and liquid samples using a custom program for the PowerViral Kit to elute 100 µL of RNA. Inhibitors were subsequently removed from the extracts using the Zymo OneStep-96 PCR Inhibitor Removal Kit (Zymo Research, Irvine, CA, USA). Eluted RNA was used immediately after extraction to quantify RNA targets of interest, and remaining RNA was stored in −80 °C. Extraction negative controls (DNA/RNA shield for solids and 1 × PBS buffer for liquids) and positive controls (BCoV spiked in DNA/RNA shield) were processed using the same protocol every time a set of samples underwent nucleic acid extraction.

Quantification

RNA targets were quantified using one-step droplet digital (dd)RT-PCR for two SARS-CoV-2 targets (N gene and S gene), BCoV, and pepper mild mottle virus (PMMoV), which acted as an endogenous internal recovery control and a fecal strength indicator (Wolfe et al., 2021a; Feng et al., 2021). A BioRad QX200 AutoDG droplet digital PCR system (Bio-Rad, Hercules, CA, USA) was used along with the BioRad One-Step RT-ddPCR Advanced Kit for Probes (Bio-Rad, Hercules, CA, USA). For SARS-CoV-2 targets, six replicate wells were used for each sample and merged in post-processing. For BCoV and PMMoV, two replicate wells were used for RNA templates diluted 1:10 in molecular grade water. Each plate included no-template PCR negative controls (water), extraction negative controls, and extraction positive controls (six replicate wells on SARS-CoV-2 plate and two on BCoV/PMMoV plate). Two positive PCR controls were run across all plates; the positive control consisted of RNA extracted from a nasopharynx swab of a high-titer COVID-19 patient from Stanford Hospital for SARS-CoV-2, direct extraction of BCoV vaccine diluted to ∼10⁶ cp/mL for BCoV, and synthetic DNA ultramer (Integrated DNA Technologies, Coralville, IA, USA) for PMMoV. Results were processed using QuantaSoft and QuantaSoft Analysis Pro (BoRad, CA) where replicate wells were merged and thresholded.

Each well needed to have at least 10,000 droplets generated to be included in the analysis. Three or more positive droplets across all replicate wells were required for a sample to be considered positive for the target. If there were less than three positive droplets, the sample was assigned as a non-detect (ND). Thresholds were chosen manually for each plate by setting a threshold for the no-template controls on the plate to have no more than two droplets above the threshold. Then the difference between this threshold and the average negative droplet fluorescence of the negative controls was set as the relative threshold difference for the plate. This difference was applied to all wells on the plate so that each well had a varying absolute threshold but a consistent relative threshold that reflected the fluctuation in the baseline negative droplet fluorescence.

The concentration per reaction was converted to copies per gram of dry weight or copies per mL of wastewater using dimensional analysis (formulas in Article S1). Errors are standard deviations as “total error” from the instrument, which includes errors associated with variability among replicate wells and the Poisson distribution. If the sample did not have a dry weight measurement associated with it due to lack of adequate solid content in the sample, an average dry weight measurement from the school’s samples was used to calculate the corresponding copies per gram of dry weight at the end of the study period.

Supplementary wastewater data

We obtained SARS-CoV-2 N gene and S gene RNA concentrations in addition to PMMoV RNA concentrations from settled solids at the wastewater treatment plant that processes the sewage from the sewershed these schools were part of. These wastewater treatment plant sample dates ranged from April 1, 2022 to June 10, 2022 to include the duration of this study. This data was acquired from the regional monitoring program to be used in a supplementary manner to examine association between school wastewater data and that of the surrounding community. These data have not been published previously but are available through the Stanford Digital Repository (https://purl.stanford.edu/km945rd8103). Methods used to acquire the regional monitoring data are described in detail by Wolfe et al (2021b) and briefly described in Article S1.

Clinical testing

Voluntary, weekly pooled clinical testing for COVID-19 was conducted at each school. Approximately 70% of the school attendees including students, teachers, and staff opted to participate in clinical testing. The schools were divided into cohorts where each cohort was tested on one day of the week (School A divided into three cohorts and School B into two). Nasal swabs of 7 to 17 individuals were collected and pooled together by the school and these samples were shipped overnight to be processed by a commercial lab using FDA-approved RT-PCR tests. The number of pools each day varied based on the number of individuals who sought testing but varied from 115 to 332 for School A and 66 to 486 for School B. Results of the pooled testing were communicated to the schools within 48 h of specimen collection.

Statistical analysis

Statistics were computed using RStudio (version 1.4.1106) and packages tidyverse, tidyr, zoo and stats. Nonparametric Mann–Whitney U test was used to test whether PMMoV RNA concentrations or PMMoV RNA concentration ratio in solid to liquid samples were significantly different between schools; associated U-statistic and p-values were reported. Fisher’s exact test was used to compare detection frequency between solid and liquid samples; associated Fisher’s exact test statistics were reported with α of 0.05 to evaluate significance of the comparison. Bonferroni correction was used only when multiple comparisons were conducted with a single hypothesis, and a new α was reported in comparisons where the correction was used.

Nonparametric Kendall’s tau was used to assess association between RNA concentrations (PMMoV and SARS-CoV-2 N, S genes) in solid versus liquid samples. To account for technical variability of wastewater measurements, a bootstrapping approach was used where Kendall’s tau was calculated using 1,000 resampling from a uniform distribution between upper and lower standard deviations of each sample concentration. Median tau and empirical p-values were calculated from this resampling; there are no other test statistic associated as this calculation involved a bootstrapping approach (Wolfe et al., 2021a). For samples that were reported as non-detects (NDs), bootstrapping bounds were defined as zero for the lower standard deviation and the lowest theoretical measurement limit for the upper standard deviation. The lowest theoretical measurement limit for each sample was calculated by determining the concentration of each sample if only three droplets across all replicate wells were positive (our criteria to be considered to have detectable concentrations of the target). Because every sample has a different number of generated droplets and dry weight associated with the sample, the lowest theoretical measurement limit varies from sample to sample (values provided in Results section).

Empirical relationships between RNA concentrations of matched solid and liquid samples were established with linear regression of log-transformed data to derive slopes and y-intercepts. Only samples that were positive for the SARS-CoV-2 target in both matched solid and liquid samples were used to investigate this relationship. Test statistics associated with linear regression such as the t-statistic along with R² value, p-value, and degrees of freedom were reported. However, it should be noted that linear regression was used mainly to derive an empirical equation for the concentrations observed, not to assess linearity of the relationship.

To compare the detection results of wastewater sampling and pooled clinical testing, we considered weekly sets of pooled clinical tests as a unit of analysis as each daily clinical pooled test corresponded to a subset of the population of the school (one cohort) but together represented the entire population that contributed clinical specimens during that week. Weekly sets of pooled clinical test results were compared to wastewater samples collected the day of or one day after specimen collection as shown in Fig. 1. This matching scheme for wastewater and clinical testing data was chosen to account for the time it takes to run clinical tests and identify an individual infected with COVID-19 (up to 48 h) to be taken out of the school system. We assumed that until identified, the infected individual was still present and contributing to the school wastewater. Therefore, detection of SARS-CoV-2 RNA in Tuesday or Wednesday wastewater samples was compared to sets of three pooled testing results at School A, done in mutually exclusive cohorts on Monday, Tuesday, and Wednesday of each week. Similarly, detection of SARS-CoV-2 RNA in Tuesday or Wednesday wastewater samples was compared to sets of two pooled testing results at School B, done in mutually exclusive cohorts on Monday and Tuesday of each week. Fisher’s exact test was used to compare solid and liquid in their ability to accurately reflect results of the positive pooled testing.

Quality assurance (QA)/quality control (QC)

Negative and positive extraction and PCR controls were ND and positive, respectively. Recovery rates of BCoV indicated no RNA extraction failures (recovery rate >10%). One sample in School A did not have a detectable concentration of PMMoV; we considered this sample to have experienced an RNA extraction failure and excluded it from further analysis. After removal of that sample, there were 91 samples (45 solid and 46 liquid samples) from the two schools remaining for inclusion in our analyses.

Measurement overview

Across solid samples, PMMoV ranged from 2.6 × 10⁶ to 1.9 ×10⁹ copies g⁻¹ dry weight (hereafter referred to as cp g⁻¹, median = 7.8 × 10⁷); across liquid samples, PMMoV ranged from 5.4 × 10² to 9.8 × 10⁵ copies mL⁻¹ wastewater (hereafter referred to as cp mL⁻¹, median = 1.4 × 10⁴) (Fig. S1). PMMoV concentrations were not different between schools for solids or liquids (two separate Mann–Whitney U tests, for solid U statistic = 181, p-value for solids = 0.10; for liquid U statistic = 281, p-value = 0.73). Median ratio of PMMoV concentrations in matched solid to liquid samples across both schools was 4.1 × 10⁴ mL g⁻¹ (n= 45, range 3.8 × 10² to 8.2 × 10⁴ mL g⁻¹). Ratios were different between schools with School B having a higher median ratio (4.8 × 10⁴ mL g⁻¹) compared to School A (3.4 × 10³ mL g⁻¹) (Mann–Whitney U test, U statistic = 156, p-value = 0.03). PMMoV concentrations in matched solid and liquid samples were positively and significantly correlated at the two schools as both aggregated data and at individual school level (Kendall’s tau > 0, empirical p-value = 0 for all, α with Bonferroni correction = 0.017; for specific values see Table 1).

In the solid and liquid fractions of each sample, N and S gene targets of SARS-CoV-2 were measured (Fig. 2). SARS-CoV-2 RNA gene concentrations in solids ranged from ND to 4.4 × 10⁴ cp g⁻¹ (N) and from ND to 1.1 × 10⁵ cp g⁻¹ (S). Across liquid samples, SARS-CoV-2 RNA concentrations ranged from ND to 41 cp mL⁻¹ (N) and from ND to 6 cp mL⁻¹ (S) (Fig. S2). The lowest theoretical measurement limit calculated for each sample ranged from 3.3 × 10³ to 1.7 × 10⁴ cp g⁻¹ (median = 9.7 × 10³) for solid samples and 0.2 cp mL⁻¹ to 0.5 cp mL⁻¹ (median = 0.3) for liquid samples (Fig. S3). Lowest observed measurements for solid samples were 8.1 × 10³ (N) and 8.9 × 10³ (S) cp g⁻¹; for liquid samples were 0.3 cp mL⁻¹ (N and S). When comparing the SARS-CoV-2 concentrations at each of the schools to that of the larger sewershed that the schools were part of only School B had a positive and significant association (Kendall’s tau = 0.23, empirical p-value = 0 for N and Kendall’s tau = 0.11, empirical p-value = 0.047 for S; Fig. S4, Table S1). All wastewater data presented in this paper can be found in the Stanford Digital Repository (https://purl.stanford.edu/sy647tw8455).

N and S were significantly correlated within both solid (R² = 0.12, slope = 0.22, p-value = 0.01, degrees of freedom = 43, t-statistic = 2.7) and liquid measurements (R² = 0.39, slope = 3.97, p-value = 2.2 × 10⁻⁶, degrees of freedom = 44, t-statistic = 5.4) (Fig. S5). Although they were positively associated, N and S measurements were not always in agreement; detection of N and S gene were in agreement (i.e., both detected or both undetected) in 62% of samples and they were in disagreement (i.e., one gene detected but the other gene undetected) in 38% of samples across both solid and liquid samples. When separated by matrix, N and S gene detection agreement remained similar (63% in liquid and 60% in solid). In subsequent analyses that involved examining the relationship between SARS-CoV-2 RNA concentration in solid and liquid samples, the two targets were analyzed independently. When comparing the detection rate of wastewater samples and clinical pooled testing, an individual solid or liquid wastewater sample was considered positive for SARS-CoV-2 if either one of the targets (S or N) was detected.

Relationship between SARS-CoV-2 RNA in solids and liquids

Out of the 45 matched liquid and solid samples, 31 solid samples were positive, and 14 samples were negative for SARS-CoV-2 RNA; 34 liquid samples were positive, and 11 samples were negative for SARS-CoV-2 RNA (Table 2). There is no significant difference between detection frequency in solid versus liquid samples (Fisher’s exact test statistic 0.64, α with Bonferroni correction = 0.017). The result did not change when samples from each school were examined separately. Solid and liquid samples agreed on detection in the majority of samples (67%) but there were samples that were positive for SARS-CoV-2 RNA in the solids but negative in the liquid and vice versa. School B samples had a higher rate of agreement (74%) between solids and liquids, but this rate was not statistically different from the rate for School A samples (59%) (Fisher’s exact test statistic 0.35).

The median ratio of SARS-CoV-2 RNA concentrations in matched solid and liquid samples where both matrices were positive for the SARS-CoV-2 N gene was 8.6 × 10³ mL g⁻¹ (n = 15, min 5.6 × 10² mL g⁻¹, max 6.5 × 10⁴ mL g⁻¹) and for the S gene was 1.6 × 10⁴ mL g⁻¹ (n = 14, min 4.0 × 10³ mL g⁻¹, max 9.3 × 10⁴ mL g⁻¹) (Table 3). The ratios of N gene and S gene are not statistically different (Mann Whitney U test, U statistic = 72, p-value = 0.16). SARS-CoV-2 RNA concentrations in matched solid and liquid samples were positively and significantly correlated for both N (median Kendall’s tau = 0.21, empirical p-value = 0) and S gene (median Kendall’s tau = 0.24, empirical p-value = 0, α with Bonferroni correction = 0.017) as aggregated data. At the individual school level, SARS-CoV-2 concentrations in matched solid and liquid samples for each school were positively and significantly correlated for School B (all median Kendall’s tau > 0, empirical p-value < 0.05, α with Bonferroni correction = 0.017; see Table 1 for exact values). When SARS-CoV-2 RNA concentration was normalized by PMMoV, SARS-CoV-2 concentrations in matched solid and liquid samples for each school were positively and significantly correlated for both targets (all median Kendall’s tau > 0, empirical p-value < 0.05, α with Bonferroni correction = 0.017; see Table 1 for exact values).

To derive an empirical relationship between the log₁₀-transformed solids and liquid concentrations, we used linear regression where y is the log₁₀-transformed solid concentration (cp g⁻¹) and x is the log₁₀-transformed liquid concentration (cp mL⁻¹), as in a Freundlich isotherm model (Schwarzenbach, Gschwend & Imboden, 2017): $C_S=K_f{C}_L^{1/n}$ where C_S is RNA concentration in the solid samples, C_L is RNA concentration in liquid samples, K_f is the Freuidlich’s constant, and 1/n is the exponent of non-linearity. For this analysis, we used only data points that had SARS-CoV-2 RNA concentration detected in both solid and liquid samples instead of substituting ND with other values. This was done to prevent the substituted values from having an impact on the derived relationship. Data from the two schools were examined in aggregate. For the N gene, the slope was 0.21 and the y-intercept was 4.1 (R² = 0.34, p-value = 0.01, degrees of freedom = 13, t-statistic for slope 2.9, t-statistic for y-intercept = 77.6), indicating n = 4.8 and K_f = 10⁴ mL g ⁻¹ in the Freundlich model. For the S gene, the slope was 0.29 and the y-intercept was 4.3 (R² = 0.11, p-value = 0.13, degrees of freedom = 12, t-statistic for slope 1.6, t-statistic for y-intercept = 57.5), indicating n = 3.4 and K_f = 10⁴ mL g ⁻¹ (Fig. 3, Table 4).

Comparison with pooled testing

Out of 16 school-weeks total (eight weeks at two schools) where pooled testing was completed, COVID-19 infection was confirmed in pooled testing for 12 of the school-weeks. During these 12 school-weeks, nine sets of the solid samples (75%) and 12 sets of the liquid samples (100%) included positive detections of SARS-CoV-2 RNA in wastewater. The difference between the detection rate of solid (75%) and liquid (100%) samples to detect SARS-CoV-2 RNA in wastewater when pooled testing was positive is not significantly different (Fisher’s exact test statistic = 0.22). In the 4 school-weeks of negative pooled clinical specimen testing, 3 weeks of solid samples (75%) and three sets of liquid samples (75%) included positive detection of SARS-CoV-2 RNA in wastewater (Fig. 4).

On Fridays, clinical specimens were not collected but wastewater samples were. Out of 15 school-days total (eight Fridays at two schools with one failed sampling event), 12 solid samples (80%) and 12 liquid samples (80%) were positive for SARS-CoV-2 RNA (Fig. S6), indicating the possibility that the week’s clinical sampling may not have removed all the individuals affected with COVID-19 from the campus.