Wastewater sequencing uncovers early, cryptic SARS-CoV-2 variant transmission

Abstract

As SARS-CoV-2 continues to spread and evolve, detecting emerging variants early is critical for public health interventions. Inferring lineage prevalence by clinical testing is infeasible at scale, especially in areas with limited resources, participation, or testing/sequencing capacity, which can also introduce biases. SARS-CoV-2 RNA concentration in wastewater successfully tracks regional infection dynamics and provides less biased abundance estimates than clinical testing. Tracking virus genomic sequences in wastewater would improve community prevalence estimates and detect emerging variants. However, two factors limit wastewater-based genomic surveillance: low-quality sequence data and inability to estimate relative lineage abundance in mixed samples. Here, we resolve these critical issues to perform a high-resolution, 295-day wastewater and clinical sequencing effort, in the controlled environment of a large university campus and the broader context of the surrounding county. We develop and deploy improved virus concentration protocols and deconvolution software that fully resolve multiple virus strains from wastewater. We detect emerging variants of concern up to 14 days earlier in wastewater samples, and identify multiple instances of virus spread not captured by clinical genomic surveillance. Our study provides a scalable solution for wastewater genomic surveillance that allows early detection of SARS-CoV-2 variants and identification of cryptic transmission.

Extended Data Figure 1:. Relationship of daily UCSD campus wastewater sampler positivity and campus clinical positives.

Black line indicates the linear fit to the data, with bootstrap 95% confidence interval shown in gray.

Extended Data Figure 2:. Relationship between genome coverage and cycle quantification values.

10x genome coverage (fraction of sites with 10 reads or greater) remains high, even for Cq values of nearly 38. Points indicate median value in each bin, while error bars indicate the median absolute deviation.

Extended Data Figure 3:. Lineage-specific prediction of variant abundance in spike-in validation samples.

A. Schematic of “spike-in” sample design. B-F. Lineage specific prediction. Proportions of each lineage in the sample are shown as a pie chart marker (Grey = Lineage A, Orange = Alpha, Pink = Beta, Turquoise = Delta, and Purple = Gamma) with error bars indicating the standard deviation from the mean, across four replicates.

Extended Data Figure 4:. Freyja more accurately estimates virus abundance, with fewer false positives.

A-B. Estimated vs expected fraction of each lineage in the mixture. The Kallisto-based approach from Baaijens et. al shows a wider range of estimates for each known mix fraction, and generally underestimates the fraction. C. False positives with abundance greater than 0.5%.

Extended Data Figure 5:. The rise of the Delta variant during Summer 2021

A. Mean SARS-CoV-2 viral gene copies/L of raw sewage (blue) collected from the Point Loma Wastewater Treatment Plant and caseload (gray) reported by the county during the same period. SARS-CoV-2 concentrations were normalized by PMMoV (pepper mild mottle virus) concentration to adjust for load changes. B. Lineage distribution in UCSD campus wastewater. C. Monthly lineage averages for wastewater collected at Point Loma Wastewater Treatment Plant during the Delta surge (N= 5, 20, 25, 7)

Extended Data Figure 6:. Quantification of deconvolution uncertainty in first detection of VOCs.

A-D. Bootstrap distributions of Freyja abundance estimates obtained by resampling read data from each sample corresponding to the first detection of that VOC in San Diego. Two samplers were found to contain Delta on the same day. First detections were also confirmed using a VOC qPCR panel, as shown in Figure 2 and Extended Data Table 3. 95% Confidence intervals for variant prevalence for each first detection event: A. Alpha: (0.232, 0.278), B. Delta: (0.336, 0.397), C. Delta: (0.676, 0.772), D. Omicron: (0.017, 0.021).

Extended Data Figure 7:. Estimated proportion of Omicron sequences in clinical data.

Omicron estimates tracked via S-gene target failure, SGTF (characteristic of Omicron lineage BA.1 and its descendants) qPCR assays for clinical samples in San Diego between November 27th, 2021-February 7th, 2022. First detection of Omicron through clinical genomic sequencing in San Diego was December 8th. Dotted line shows a rolling average with a window size of seven days.

Figure 1:. Campus sampling locations and SARS-CoV-2 testing statistics.

A. Geospatial distribution of the 131 actively deployed wastewater autosamplers and the corresponding 360 university buildings on the campus sewer network. Building-specific data have been de-identified in accordance with university reporting policies. B. Campus wastewater and diagnostic testing statistics over the 295 day sampling period (WW = wastewater, positivity is the fraction of WW samplers with a positive qPCR signal). C.Virus diversity in wastewater and clinical samples: Boxplots of Shannon entropy (top) and richness (bottom) for each sample type.

Figure 2:. Sample deconvolution robustly recovers relative virus abundance.

A. Subset of lineage defining mutation “barcode” matrix. Each row represents one lineage (out of >1000 lineages included in the UShER global phylogenetic tree), and individual nucleotide mutations are represented as columns. B. Single nucleotide variant frequencies obtained from iVar used for recovering relative abundance of each lineage. C. Schematic of the spike-in validation experiment. D. Depth-weighted de-mixing estimates of the virus abundance versus expected/known abundance. Details on lineage specific predictions are provided in Extended Data Figure 3. E. Comparison of wastewater sample deconvolution with VOC qPCR panel, with lookup table (bottom) showing amino acid mutations corresponding to each variant.

Figure 3:. Freyja recovers early and cryptic transmission of SARS-CoV-2 variants of concern

A. Timeline and normalized epidemiological curves for VOC detection in both wastewater and clinical sequences from San Diego County for the 3 major VOCs in circulation during the sampling period. Both Alpha and Delta are detected first in wastewater before clinical samples. Markers for clinical detections correspond to the ceiling of the daily detection count divided by 30 (e.g. 1–30 samples= one marker, 31–60 = two markers), while wastewater markers correspond to a single detection. B. Timeline and epidemiological curves for VOC detection in the campus samples. Markers correspond to a single detection event for both clinical and wastewater surveillance. All wastewater detections correspond to an estimated VOC prevalence of at least 10%.

Figure 4:. Deconvolution recovers a fine-grained estimate of virus population dynamics.

A. Prevalence of SARS-CoV-2 variants in UCSD clinical surveillance, and B. Variant prevalence in all clinical samples collected in San Diego County. C,D. Variant prevalence in wastewater at UCSD as well as the greater San Diego County (includes wastewater samples collected from Point Loma wastewater treatment plant as well as public schools in the San Diego districts). Further analysis of Point Loma wastewater samples is shown in Extended Data Figure 5. All curves show rolling average, window ±10 days. “Other” contains all lineages not designated as VOCs. Bottom panels show number of sequenced samples per day.

Figure 5:. Community wastewater enables early Omicron detection and reveals lineage dynamics.

A. Prevalence of SARS-CoV-2 VOCs in wastewater collected from the Point Loma wastewater treatment plant from late September 2021 to early February 2022. B. Estimated VOC concentrations, prevalence estimates scaled by normalized viral load in wastewater. C,D. Lineage-specific estimates of prevalence and concentration. All curves show an adaptive rolling average calculated using a local linear approximation (Savitzky-Golay filter) of virus copies/L, with window size ± 1 sampling date.

Figure 6:. Wastewater identifies clinically known and unknown virus transmission.

A-C. Maximum likelihood phylogenetic trees for each of the dominant variants of concern using high quality samples obtained at UCSD, as well as a representative set of sequences from the entire United States. Wastewater sequences from the same sampler that differ by 1 or fewer SNPs are denoted with a red asterisk. For all sequences, consensus bases were called at sites with >50% nucleotide frequency. Location information is provided for select outbreaks. D. Pairwise comparison of collection date for matching and near-matching wastewater and nasal swab samples obtained at UCSD. Positive values indicate earlier collection in nasal swabs, and negative values indicate earlier detection in wastewater.