Wastewater sampling

Wastewater samples were collected from a WWTP located in Suffolk County, NY. The WWTP serves ~330,000 people in its sewer catchment and treats ~30.5 million gallons of wastewater daily. Untreated raw sewage influent was collected via autosampler every 15 min to make up a 24-h composite sample, refrigerated upon collection to <6 °C. The 24-h composite sample was subsampled into a 500-ml polypropylene bottle, stored in a cooler with ice packs, and then transferred to our lab (~1 h drive) for subsequent analyses. Viral analysis was performed immediately upon sample receipt and the entire procedure was completed within a day. The remaining sample was then split into two aliquots and stored at −80 °C without adding any preservatives, one for drug analysis and the other as an archived sample. Suspended particles in the samples for drug analysis were removed through vacuum filtration (1 µm glass fiber) prior to freezing. The subsequent drug analysis was performed within 3 weeks. Sampling was initiated on June 2020, with daily sampling from June 3 to June 9, weekly sampling from June 9 to July 7, and biweekly sampling from July 7 to December 22, 2020, and twice weekly sampling from January 2021 through January 6, 2022.

Detection and quantification of SARS-CoV-2 RNA

Twenty-four-hour composite samples of raw sewage were centrifuged at 4200 rpm for 30 min at 4 °C in order to remove large particles and debris before polyethylene glycol (PEG) precipitation. To evaluate the viral recovery rates from wastewater, bovine coronavirus (BCoV), which belongs to the same genus as SARS-CoV-2, was spiked into the supernatant. The viral particles in 40 ml of samples were precipitated with PEG 8000 (Millipore Sigma, Burlington, MA) and NaCl (5 M, Millipore Sigma, Burlington, MA) and then incubated overnight at 4 °C. RNA from the PEG-precipitated wastewater was extracted by Qiagen QIAamp DSP viral RNA mini kit (Qiagen, Hilden, Germany) according to manufacturer’s instructions and eluted in 100 µl by nuclease-free water. The concentrations of RNA were measured by NanoDrop One Spectrophotometer (Thermo Fisher Scientific, Waltham, MA). All RNA samples were stored at −80 °C and subjected to cDNA synthesis within the same day of RNA extraction to avoid losses associated with storing and freezing and thawing RNA extracts.

Reverse transcription was performed by High Capacity RNA-to-cDNA Kit (Applied Biosystems, Waltham, MA) at 37 °C for 60 min, and stored at −20 °C until further analysis. The cycling condition was 95 °C for 10 min, followed by 40 cycles of 95 °C for 5 s and 55 °C for 40 s, and 98 °C for 10 min. The total volume of each reaction was 14.5 µl containing 7.25 µl of QuantStudio 3D Digital PCR Master mix v2 (Applied Biosystems, Massachusetts, USA), 0.725 µl of primer and probe (N1/ BCoV), 0.725 µl of TaqMan® Copy Number Reference Assay RNase P (as an internal control, Applied Biosystems, Waltham, MA), 4.8 µl of nuclease-free water, and 1 µl of cDNA template. Digital PCR was performed using N1 primers and probe set from 2019-nCoV RUO Kit (IDT # 10006713) with the CDC-recommended sequence and BCoV set against the BCoV gene as an external reference on a QuantStudio 3D Digital PCR (Applied Biosystems, Massachusetts, USA). Nuclease-free water was used as non-template control (NTC) and plasmids containing the complete nucleocapsid gene from 2019-nCoV (IDT # 10006625) were used as a positive control. Data analysis was performed with the online version of the QuantStudio 3D AnalysisSuite Cloud Software. The limit of detection for the N1 gene was 1.4 copies/reaction.

Detection of COVID-19 treatment drugs and other pharmaceuticals

Standards of COVID-19 treatment drugs (i.e., remdesivir, chloroquine, and hydroxychloroquine), other pharmaceuticals (OTC drugs), and their corresponding isotopically labeled compounds were purchased from Toronto Research Chemicals Inc (Ontario, Canada), Sigma-Aldrich (MO, USA), Fisher Scientific (MA, USA), Cerilliant (TX, USA), CDN Isotopes (Quebec, Canada), and Cambridge Isotope Laboratories (MA, USA). A list of all compounds used in this study is presented in Table S1. The 26 pharmaceuticals, including OTC drugs and the metabolites, were chosen because of their relatively high prescriptions per population in the U.S. and high environmental detection frequencies in previous studies [38, 39].

Due to the low concentration of COVID-19 treatment drugs, SPE was required to concentrate the sample for detection. In contrast, other pharmaceuticals were measured by direct injection after dilution. In brief, a 100-ml sample was transferred out and spiked with a surrogate standard (hydroxychloroquine-D4) to trace the extraction yield prior to SPE. After conditioning the SPE cartridge (Waters Oasis HLB, 200 mg, 6 cc) with methanol and deionized water, the whole sample was loaded onto the cartridge, and after which the cartridge was eluted sequentially with 4 ml methanol, resulting in ~25-fold preconcentration. Extracts were stored at −20 °C until analysis. Prior to analysis, the extract was diluted with deionized water (50:50 MeOH: H2O) and spiked with the internal standard. For other pharmaceuticals, a 100-µl sample was taken out and diluted 10-times with deionized water and methanol to constitute a final concentration of 10% methanol. The isotopically labeled internal standards were then added before analysis. The detailed information for the surrogates and internal standards is shown in Table S1.

Detection and quantification of the target compounds in extracts were carried out using an Agilent 6495B triple-quadrupole mass spectrometer (LC-MS/MS) with an electrospray ionization source in positive ion mode (ESI+), using Multiple Reaction Monitoring (MRM) to monitor the precursor ions and product ions (Table S1). The detailed instrumental conditions are shown in Table S2.

Stability of COVID-19 treatment and OTC drugs in wastewater

The wastewater temperature in underground sewer pipes in the study area was measured to vary between 10–12 °C in winter and 18–20 °C in summer, as provided by wastewater operators. The travel time of sewage from houses to the WWTP in the study area ranged from 40 min to 8 h, depending on the distance. Additionally, the collected wastewater sample could reside in the composite sampler for up to 24 h at 6 °C. Thus, we performed a controlled experiment to assess the stability of the analytes of interest at different temperatures within 24 h. In brief, a suite of COVID-19 treatment and OTC drugs was spiked (50 ng each) into 50 ml of unfiltered raw wastewater. The spiked wastewater samples were stored at 4, 12, and 20 °C. Each temperature had triplicate samples. At t = 0 h and t = 24 h, 5 ml aliquot was taken out for analysis and followed the sample preparation procedure described above for LC-MS/MS analysis.

COVID-19 cases

New reported confirmed cases of COVID-19 were recorded by the Suffolk County Department of Health, NY. Data at the zip code level were shared with our research team to support COVID-19 research in the region. We identified the 13 zip codes in the catchment area of the WWTP and summed the number of cases daily in the catchment area to create a 7-day rolling average number of cases.

We also received reports from Stony Brook University Hospital of daily hospitalized cases, and milligram of COVID-19 treatment drugs, hydroxychloroquine and remdesivir prescribed daily beginning Oct 3, 2020 to the present. This hospital is not physically located in the catchment area but is the closest level 1 Trauma center to the catchment area, and receives patients from the catchment area. These data, therefore, are not used as proxies for the amount of remdesivir or hydroxychloroquine in the catchment area, but rather are useful for understanding temporal trends in prescriptions of these treatment drugs in the region. As shown in Fig. S1, remdesivir usage can reflect the case trend in the hospital, whereas hydroxychloroquine usage remains relatively stable over time.

Population correction

Estimating the actual population contributing to sewage flow during the sampling period is challenging but is essential as it directly influences the concentration of biomarkers in wastewater. Several endogenous and exogenous human biomarkers have been proposed to serve as a tool for population normalization [16, 40]. In this study, we selected caffeine, a stimulant excreted in human urine, to estimate the serviced population in the WWTP sewershed because its level in wastewater is known to be stable and features <10% degradation within 24 h as shown in our preliminary experiment (Fig. S2). Over the ~20-month period, the caffeine concentration in the samples showed little variation over time with a mean concentration of 88.2 ± 20 µg/l (range: 48.2–148 µg/l). Time series data of SARS-CoV-2 RNA, COVID-19 treatment drugs, and other pharmaceuticals were normalized by caffeine using the equation below:

$$\frac}}}}}}\,}}}}}}\,}}}}}}]}_}}}}}}}]}_}}}}}}}]}_}}}}}}}}}$$

(1)

where [Virus or drug]t is the virus or drug concentration at time = t, [Caffeine]t is the caffeine concentration at time = t, and [Caffeine]avg is the average caffeine concentration.

Model development

We developed several Bayesian models to predict confirmed cases with concentrations of virus gene copies and other biomarkers in wastewater samples. The Bayesian framework, as compared to classical statistics, allows us to update our current models with future data collection. Except for the confirmed cases, all other variables were adjusted in the following ways for ease of modeling or interpretation. First, virus concentration was log-transformed with base 10. All other measured chemicals were divided by their maximum value in the sample and rescaled into values between 0 and 1. This way, we retained zero as a reference point while being able to make sense of the priors across variables in our modeling. In order to select our variables of interest to fit the model, data exploration was performed and described in the Supplementary Information.

We regressed our dependent variable, the confirmed cases at time t (Ct), on lags of the predictor variable(s) (X1, t−1, X1, t−2, X1, t−3,…; X2, t−1, X2, t−2, X2, t−3,…; X3, t−1, X3, t–2, X3, t−3,…, etc.,), via a generalized linear model (GLM). The number of confirmed cases was modeled as a binomial distribution with the trial number equal to the population in the sewershed area (N = 3.3 × 105). The probability of individual infection is a logistic function of linear combinations of our predictor variables. The models have a general form:

$$_ \sim }}}}}}\left(N,p\right)$$

(2)

$$}}}}}}\left(\right)=\alpha +\mathop\limits_=1}^\left(\mathop\limits_=1}^_}_\right)$$

(3)

$$\alpha \sim }}}}}}\left(\bar,_\right)$$

(4)

$$_} \sim }}}}}}\left(}_},__}}\right)$$

(5)

where m = number of substances (e.g., viral or/and chemical concentrations) and n = number of lags (observations before the focal day) used for prediction. We considered three specific sets of models. For each set of models, we tested possible combinations of variables according to Data Exploration and Cross Correlation (see SI). The choice of priors was examined by predictive simulations in the SI, and the posterior distributions of parameters were estimated using Markov Chain Monte Carlo (MCMC). We retained the model with the best predictive performance based on the Watanabe–Akaike Information Criterion (WAIC). Two general rules were also applied in the modeling. First, we used consecutive lags because the change of predictor variables was more likely to have gradual effects on our outcome variable. Second, to avoid overfitting, we maintained at least ten observations for each predictor included in the models, as 111 observations were present in our sample.