Read-level mutation linkage in wastewater sequencing reveals cryptic virus evolution and persistence

Read-level mutation linkage in wastewater sequencing reveals cryptic virus evolution and persistence Dylan A. Pilz1, Christine M. Aceves1, Leandro Collares1, Karthik Gangavarapu1, Joshua I. Levy1#*, Kristian G. Andersen1# 1. Department of Translational Medicine, The Scripps Research Institute, La Jolla, CA, USA # Senior author * Corresponding author Email: [email protected] (JIL) Abstract Wastewater sequencing has become widely adopted as a complement to clinical genomic surveillance, providing population-level monitoring capable of detecting “cryptic variants”, viruses that are rarely, if ever, observed clinically. However, current wastewater analyses are limited by a lack of integration with up-to-date surveillance data and systematic frameworks for tracking virus evolution. Here, we develop an end-to-end workflow for longitudinal cryptic variant detection and tracking that leverages global SARS-CoV-2 genomic surveillance. Applying this approach for a two-year period, we identified more than 2,000 cryptic variants and characterized their emergence and persistence. Using graph analysis, we inferred relationships between cryptic variants and their likely ancestors, uncovering multistep evolutionary trajectories. Our workflow provides a scalable, accessible, and real-time solution for monitoring cryptic variants and community-level virus diversity. Introduction During an outbreak, timely monitoring of virus spread and evolution is critical to inform public health guidance and interventions1,2. Clinical genomic surveillance can effectively track transmission but is resource-intensive and biased by health system access and reporting disparities3,4. To strengthen surveillance efforts, wastewater has become widely adopted for population-level monitoring with less susceptibility to sampling bias5–8. Wastewater surveillance closely tracks clinical estimates of disease burden and variant prevalence and often enables variant detection prior to clinical surveillance9–11. However, current analysis frameworks have only begun to characterize the extent of pathogen diversity and evolution observable in wastewater, much of which is rarely captured through clinical surveillance12,13. Because short read metagenomic sequencing of complex intraspecies mixtures including wastewater often leads to ambiguous sequence reconstruction14, wastewater-specific methods leveraging prior lineage annotation have become widely used for tracking lineage prevalence15–18. While enabling robust tracking of lineages, there remains a need for wastewater tools that can detect, characterize, and track virus spread and evolution, especially for emerging viruses. At the extreme, some viruses have been observed in wastewater for months without detection via clinical surveillance19,20. These “cryptic variants” have been implicated in the emergence of many SARS-CoV-2 variants of concern20,21, including through accelerated accumulation of virus mutations during chronic infections22–24. However, current computational tools for the detection and tracking of cryptic variants12,25 are limited in their computational efficiency, access to up-to-date surveillance data, and scalability. Enabling real-time, comprehensive cryptic variant tracking will require bioinformatic tools that are independent of cladistic nomenclature and integrated with global genomic surveillance. Here, we show that a workflow for longitudinal monitoring of virus diversity using mutation-linkage in wastewater sequencing, integrated with global genomic surveillance, enables tracking of cryptic variant evolution and spread. To recover read-level linkage information from wastewater and assess its significance, we developed a bioinformatic tool, coVar, to identify linked mutations and used outbreak.info26 to query known clinical sequence diversity. Deploying this approach in San Diego over a two-year period, we found that the frequency of cryptic variants followed increases in virus abundance during the early stages of lineage waves. We observed rapid virus diversification following lineage introductions, with some cryptic descendants persisting for months undetected by clinical surveillance. Using shortest-path graph analysis, we recovered stepwise evolutionary histories of circulating variants and tracked their geographic and temporal spread. Our approach provides an efficient, real-time, scalable, and systematic solution for tracking cryptic variant evolution and spread. Results Contextualizing mutation-linkage with real-time global genomic surveillance While clinical surveillance relies on single-timepoint sampling of a subset of symptomatic individuals, wastewater sampling can capture all virus shedding in the community independent of reporting and patient symptoms (Fig. 1A). To characterize virus diversity via wastewater sequencing, we developed coVar, a bioinformatic tool for quantifying the prevalence of linked mutations in virus sequencing of mixed samples. coVar recovers co-occurring mutations on sequencing read pairs in a specified genomic region, with a 5.26x decrease in CPU time over available tools12,25 in a single threaded test (Supp. Fig. 1A), further improved to 12.57x with multi-threading in a simulated sample (Supp. Fig. 1B), while closely matching results of leading tools12,25 (Supp. Fig. 1C). Figure 1: Clinical-integrated wastewater monitoring workflows enable real-time tracking of community virus spread and evolution A. Clinical surveillance includes a subset of all cases, sequences a fraction of those cases, and relies primarily on single-timepoint testing, preventing longitudinal characterization of chronic infections. Wastewater surveillance enables aggregate longitudinal measurement of all viruses circulating in the community. B. Types of viruses present in wastewater sequencing, including designated lineages, known sub-lineages spreading in the community, and clinically unobserved “cryptic” pathogen lineages. C. Steps involved in linked mutation cluster monitoring and detection of local diversity and cryptic sequences. Sequencing reads are aligned to a reference, after which Freyja is used to tabulate linked mutation clusters in the sample. The outbreak.info mutation prevalence API is used to determine how frequently a particular cluster of mutations has been seen in clinical sequences. Contextualizing read-linked mutations with global genomic surveillance can inform our understanding of community virus circulation, but in practice, attributing sequencing reads to specific lineages is challenging, particularly for short read sequencing (Fig. 1B). To identify potential virus lineages associated with linked mutations and distinguish between known and cryptic diversity, we built a workflow using the outbreak.info API26 to query global clinical sequence data for viruses possessing the identified mutation-linkage (Fig. 1C). Overall, this end-to-end workflow dramatically reduces resource and expertise requirements needed to perform these analyses and easily scales to pandemic-scale sequencing data without requiring users to process it themselves. The workflow is available as a containerized nextflow pipeline and can easily be run on a personal laptop. Wastewater reveals the timing and ancestors of cryptic variants Using this workflow, we examined the spike (S) gene of 946 samples collected from 3 San Diego wastewater treatment plants during January 2023-December 2024 (Supp. Fig. 2A,B). Among unique linked-mutations observed in San Diego wastewater, we found 2,304 that were present in <10 sequences globally, which we defined as “cryptic variants”, including 932 that were not detected at all. Comparing across samples, we found that cryptic variant detection depended on sequencing depth (Supp. Fig. 3). To control for potential confounds including false positive detection due to differences in sequencing depth, we normalized variant counts by read depth in subsequent analyses. To investigate the determinants of cryptic variant detection, we first considered virus lineage dynamics during the study period (Fig. 2A). We found that cryptic variant detections, normalized by sequencing depth, increased following major lineage introductions and during lineage sweeps, often for weeks to months after viral loads started to decrease (Fig. 2A,B). For example, for the JN.1 lineage, which drove the highest observed viral loads over the study period, we observed a 5.5x increase in cryptic variant detections from November 2023 until January 2024, which was sustained for 6 weeks. Between JN.1 reaching 20% prevalence on December 3rd, 2023 and reaching 90% prevalence 6 weeks later, we observed 34% of all normalized cryptic variant detections. To attribute cryptic variants to specific parent lineages, we identified all lineages containing a subset of the observed mutations based on existing lineage definitions9 (Fig. 2C). In cases with multiple potential parent lineages, we assigned variants to the most recent common ancestor (MRCA). Using this approach, we found that a small subset of cryptic variants (0.7%) could be attributed to specific parent lineages, including XBB.1.5, XBB.1.16, XBB.1.9, JN.1, and XEC. The remaining cryptic variants were only identifiable as descendants of basal lineages (e.g., lineage B.1). We found that lineage-attributable cryptic variants could often be detected shortly after lineage detection via wastewater. For many lineages, including XBB.1.16, JN.1, and XEC, we observed multiple cryptic descendants prior to the parent lineage reaching 5% prevalence (Fig 2A,C). However, we detected lineage-attributable variants more consistently once their parent lineages became dominant. These detections were remarkably persistent, with XBB.1.5 and JN.1 cryptic descendants each detected for more than 6 months after their parent became the dominant lineage (Fig. 2C). Figure 2: Cryptic variants emerge following lineage introduction and diversification A. Relative SARS-CoV-2 lineage prevalences observed in Point Loma wastewater treatment plant (WWTP), processed using Freyja. B. Log of normalized cryptic counts per month (blue). Briefly, cryptic variant occurrences are summed by month across all three WWTPs, before dividing by the median per-site sequencing depth of all samples collected that week. The log average viral load (genome copies / L) is shown in black. Light and dark gray bars showing periods of 5% and 20% VOC prevalence, respectively. C. Normalized cryptic counts broken down by lineage using Freyja barcodes. Extended persistence of cryptic descendants following lineage introduction To track cryptic variants and their persistence, we quantified both the number of detections in wastewater and the time between the first and last detection. Excluding cryptic variants that were only observed once, we found that 72% of cryptic variants had 10 or fewer occurrences, and 87% had 50 or fewer occurrences (Fig. 3A). 75% were detected for more than 100 days, and 2% were detected for more than 500 days (Fig. 3B). Figure 3: Mutation-linkage network analysis identifies putative evolutionary steps underlying cryptic variants, and tracks the spectrum of community virus diversity A. Histogram showing distribution of cryptic variant detection counts. B. Histogram showing the time between the first and last detection for each cryptic variant. The cumulative distributions for each are shown on the right axes. C. Graph showing the rapid diversification starting from XBB.1.5-derived lineage, where nodes represent observed variants, and edges indicate acquired mutations. The total number of wastewater detections is shown in white, the number of global clinical detections is in red, and box shading indicates WWTPs where a variant was observed. For each of the 10 cryptic descendants shown, the sampling location of individual detections is colored by WWTP. D. Examples of cryptic, stepwise evolution of variants. Dates for which a given variant was first and last detected are shown above each node. To track the evolution and geographic distribution of cryptic descendants, we used a graph-based approach to relate observed linked mutations to likely parent lineages. For XBB.1.5, for example, we consistently detected cryptic descendants throughout 2023, including some restricted to an individual WWTP, while others were found in all three catchments (Fig. 3C). The most persistent of these cryptic variants, which acquired S:K462L and S:DEL463, was detected 29 times over a 12-month window and was observed at all three WWTPs. Applying this approach to other lineage introductions, we found that the number of unique descendants increased roughly linearly over time, with short periods of accelerated detection Graph analysis recovers the putative stepwise evolution of cryptic variants To understand the evolutionary sequence preceding cryptic variant emergence, we used a graph-based shortest-path method to infer evolutionary paths through likely intermediate sequences observed in wastewater. In many instances, cryptic variants acquired several mutations over the course of several distinct evolutionary steps. For example, prior to the cryptic variant with S:T19F, S:DEL25/27, S:DEL69/70, S:S50L, S:R21T, S:DEL31, we identified likely co-circulation of S:T19F, S:DEL25/27, and S:DEL69/70 with its descendants that acquired S:S50L, and others that also gained S:R21T. In January 2024, we observed that one of these descendants with S:R21T had also acquired S:T20V, and then S:DEL31 in May 2024. (Fig. 3D, top). We observed many virus sequences and likely descendants in the same samples, indicating that viruses often co-circulated at low frequency with their descendants. In many cases, the parent lineage persisted longer than its descendants, indicating that the acquired mutations were not sufficiently advantageous for the descendant to outcompete its parent or other circulating lineages. For instance, the mutation-linkage S:N679K, S:S711P, S:P681H was detected from January 2023 until January 2024, however, and acquired S:A701V in late August 2023 (Fig. 3D, middle). However, this descendant was only detected until November 2023. Similarly, mutation-linkage S:P139Q, S:G142D was consistently detected from January to December of 2023, but its cryptic descendant with S:H146Q was only observed until October 2023 (Fig. 3D, bottom). Discussion By combining read-level mutation-linkage analysis of wastewater sequencing using coVar with real-time pandemic-scale genomic surveillance with outbreak.info, we developed an efficient and scalable method to detect and monitor cryptic variant evolution, spread, and persistence. We found that cryptic variants typically emerged shortly after lineage introductions and became more frequent after increases in wastewater viral loads. These cryptic variants could often be detected for months, with 54% detected for more than 100 days. Using a graph-based inference approach, we recovered step-wise evolutionary trajectories among these persistent cryptic descendants, and tracked their circulation and that of their parents over time. Due to the complex intra-species mixtures found in wastewater, achieving the full potential of wastewater genomic surveillance is likely to require a combination of approaches. Tools including Freyja15 and Lollipop27 are effective for inferring lineage prevalence from wastewater, but have limited utility for the identification of novel lineages. The recently proposed WEPP tool aims to leverage phylogenetic placement for haplotype recovery from mixtures, but it is not suitable for identifying low-frequency variants 28. Our approach builds from earlier methods for detecting cryptic variants12,25, providing an accessible and computationally efficient workflow for characterizing virus diversity and evolution using mutation-linkage. Integration of our coVar method with outbreak.info enables streamlined querying of global genomic surveillance without requiring custom lineage definitions, database infrastructure, or complex setup. While the observed persistence and stepwise accumulation of mutations in wastewater aligns with published findings on chronic infections in immunocompromised hosts, these chronic cases represent only a small fraction of cases. Biases in clinical sampling, however, which can miss localized transmission, asymptomatic infections, and communities with limited access to testing, are likely the primary driver of wastewater-only virus detection. Accordingly, wastewater may be a valuable tool to better understand these structural biases and guide public health resource allocation. Wastewater offers a critical population-level complement to case-based genomic surveillance but also presents challenges including limited spatial resolution, fragmented genomic sequences, and variable detection sensitivity. While these factors can complicate functional interpretation for public health responses, integrating wastewater findings with clinical sequencing, case reporting, and serological surveys has the potential to support more comprehensive and contextually aware surveillance. Technical modifications including deeper sequencing can further improve sensitivity for low-frequency circulating cryptic variants, albeit generally with increased cost and higher false positive risk. Together, our workflow lowers the barrier to characterizing community virus diversity in wastewater and enables detailed exploration of the timing, evolution, and persistence of cryptic variants. It enables real-time, efficient tracking of circulating viruses, expanding the utility of wastewater sequencing as a public health tool. Methods Physically-linked mutation extraction and cryptic variant detection In order to detect physically linked mutations, coVar (https://github.com/andersen-lab/covar) iterates over all reads in a provided alignment file using rust-htslib. The tool then identifies single nucleotide variants (SNVs), insertions, and deletions present on the same read pair, and outputs frequency and coverage information of each unique linked mutation cluster. coVar was run using minimum base quality 20 and minimum sequencing depth 10, over the SARS-CoV-2 spike gene. Mutation clusters detected in at least 2 distinct samples, having at least 2 co-occurring mutations, were kept. To detect cryptic mutation clusters, clinical prevalence of each cluster was ascertained via the outbreak.info API, and clusters with fewer than 10 clinical detections worldwide over the SARS-CoV-2 pandemic were considered cryptic. This workflow is also available as a nextflow pipeline (https://github.com/dylanpilz/cryVar). Comparison between variant prevalence, cryptic variant abundance, and viral load SARS-CoV-2 variant prevalences are resolved from wastewater using Freyja (https://github.com/andersen-lab/Freyja) and smoothed using a non-uniform Savitzky-Golay filter. To account for differences in sequencing depth when counting cryptic detections, we bin the cryptic counts by week before normalizing by the median sequencing depth for the samples within that week. Viral load is determined by averaging the mean viral load (genome copies/L), by week, from Encina, Point Loma, and South Bay WWTPs. Both cryptic counts and viral loads are smoothed using a 4-week rolling average filter. Cryptic variant clusters were assigned to candidate parent lineages by matching linked SNVs to the Freyja SARS-CoV-2 barcodes, using the longest common substring of the de-aliased PANGO lineage in order to handle ties. Note that in most instances, it was not possible to recover the parent lineage due to a lack of informative SNVs on many cryptic mutation clusters. Cryptic variant persistence and diversification The number of cryptic detections as well as the span of days between the first and last detection of each cryptic variant were summarized in histograms using 30 bins. Cumulative distributions are overlaid on each histogram, normalized to 1. Network analyses of cryptic evolution To identify diversification of cryptic mutation clusters from non-cryptic parent clusters, parent clusters are manually curated, before identifying cryptic supersets of each parent cluster. The longest branch length is limited to 1, and the maximum outdegree is limited to 10 for simplicity. In order to determine the stepwise accumulation of mutations leading to cryptic variants, mutational precursors are determined backwards by identifying subsets of mutations from a cryptic mutation cluster. These parent-child relationships are then constructed into a DAG, where nodes represent observed haplotypes, and edges represent acquired mutations. All data, analyses, as well as figure generation scripts can be found at https://github.com/andersen-lab/covariants-paper Acknowledgements We would like to thank all of the members of Modjadji wastewater surveillance initiative and the San Diego Epidemiology and Research for COVID Health (SEARCH) alliance. This work has been funded by Centers for Disease Control and Prevention (75D30120C09795 KGA), National Institutes of Health (5T32AI007244-38 to JIL; 3U19AI135995-03S2 to KGA; U19AI135995 to KGA; U01AI151812 to KGA, UL1TR002550 to KGA), the Gates Foundation (057213 to JIL, KGA), and the Conrad Prebys Foundation (to KGA). Ethics declarations KGA has received consulting fees for advising on SARS-CoV-2, variants and the COVID-19 pandemic. The other authors declare no competing interests. References 1. Kamalrathne, T., Amaratunga, D., Haigh, R. & Kodituwakku, L. Need for effective