Genetic Diversity and Evolutionary Convergence of Cryptic SARS-CoV-2 Lineages Detected Via Wastewater Sequencing.

Wastewater-based epidemiology (WBE) is an effective way of tracking the appearance and spread of SARS-COV-2 lineages through communities. Beginning in early 2021, we implemented a targeted approach to amplify and sequence the receptor binding domain (RBD) of SARS-COV-2 to characterize viral lineages present in sewersheds. Over the course of 2021, we reproducibly detected multiple SARS-COV-2 RBD lineages that have never been observed in patient samples in 9 sewersheds located in 3 states in the USA. These cryptic lineages contained between 4 to 24 amino acid substitutions in the RBD and were observed intermittently in the sewersheds in which they were found for as long as 14 months. Many of the amino acid substitutions in these lineages occurred at residues also mutated in the Omicron variant of concern (VOC), often with the same substitution. One of the sewersheds contained a lineage that appeared to be derived from the Alpha VOC, but the majority of the lineages appeared to be derived from pre-VOC SARS-COV-2 lineages. Specifically, several of the cryptic lineages from New York City appeared to be derived from a common ancestor that most likely diverged in early 2020. While the source of these cryptic lineages has not been resolved, it seems increasingly likely that they were derived from immunocompromised patients or animal reservoirs. Our findings demonstrate that SARS-COV-2 genetic diversity is greater than what is commonly observed through routine SARS-CoV-2 surveillance. Wastewater sampling may more fully capture SARS-CoV-2 genetic diversity than patient sampling and could reveal new VOCs before they emerge in the wider human population. Author Summary: During the COVID-19 pandemic, wastewater-based epidemiology has become an effective public health tool. Because many infected individuals shed SARS-CoV-2 in feces, wastewater has been monitored to reveal infection trends in the sewersheds from which the samples were derived. Here we report novel SARS-CoV-2 lineages in wastewater samples obtained from 3 different states in the USA. These lineages appeared in specific sewersheds intermittently over periods of up to 14 months, but generally have not been detected beyond the sewersheds in which they were initially found. Many of these lineages may have diverged in early 2020. Although these lineages share considerable overlap with each other, they have never been observed in patients anywhere in the world. While the wastewater lineages have similarities with lineages observed in long-term infections of immunocompromised patients, animal reservoirs cannot be ruled out as a potential source.

Introduction

SARS-CoV-2 is shed in feces of infected individuals [1, 2], and SARS-CoV-2 RNA can be extracted and quantified from community wastewater to provide estimates of SARS-CoV-2 community prevalence [3, 4]. This approach is especially powerful since it randomly samples all community members and can detect viruses shed by individuals whose infections are not recorded, such as asymptomatic individuals, those who abstain from testing, or those who test at home [5, 6]. Additionally, SARS-CoV-2 RNA isolated from wastewater can be sequenced using high-throughput sequencing technologies to define the composition of variants in the community [7-9].

The continuing evolution of SARS-CoV-2 [10] and the appearance of variants of concern (VOC), such as the Omicron VOC [11], highlight the importance of maintaining a vigilant watch for the emergence of unexpected, novel variants. The fact that the origins and early spread of the Alpha and Omicron VOCs were not observed strongly motivates efforts to detect and monitor novel variants [12]. However, whole genome sequencing of SARS-CoV-2 RNA isolated from wastewater often suffers from low sequencing depth of coverage in epidemiologically relevant areas of the genome, such as the Spike receptor binding domain (RBD)[13-15]. Additionally, because wastewater may contain a mixture of viral lineages and whole genome sequencing relies on sequencing small fragments of the genome, computational strategies to identify variants with linked mutations often fail to identify lineages present at low concentrations [16]. These features have made it difficult to detect unexpected, novel variants from wastewater samples from whole genome sequencing data.

To address these issues, we developed a “targeted” sequencing approach that amplifies and sequences the Spike RBD of the SARS-CoV-2 genome as a single amplicon (Fig. 1A) [8, 9]. Since the Spike RBD is relevant to SARS-CoV-2 infectivity, transmission, and antibody-mediated neutralization [17-21], this approach ensures that the RBD receives high sequencing coverage. Additionally, RBD sequencing enables linkage of polymorphisms, forming short, phased haplotypes [16]. These phased haplotypes permit easier lineage identification, even at low concentrations, if the targeted sequence(s) are rich in lineage-defining polymorphisms [9].

Fig. 1

RBD amplification. A. Schematic of regions targeted by the RBD and S1 primer sets (see Methods for primer sequences). Overview of the SARS-COV-2 Spike RBD lineages identified in B. the MO33 sewershed and C. the MO45 sewershed. Each row represents a unique lineage and each column is an amino acid position in the Spike protein (left). Amino acid changes similar to (green boxes) or identical to (orange boxes) changes in Omicron (BA.1) are indicated. The heatmap (right) illustrates lineage (row) detection by date (column), colored by the log10 percent relative abundance of that lineage.

Using our targeted sequencing approach, we identified and previously reported circulating VOCs in different sewersheds around the United States [8, 9]. Variant frequencies in these sewersheds closely tracked VOCs frequency estimates from clinical sampling in the same areas [8, 9]. However, in some locations, we noted the presence of cryptic lineages not observed in clinical samples anywhere in the world. Several of these lineages contained amino acid substitutions rarely reported in global databases such as gisaid.org [22-24] (e.g., N460K, Q493K, Q498Y, and N501S) [8]. Interestingly, polymorphisms in these lineages show considerable overlap with the Omicron VOC, suggesting convergent evolution due to similar selective pressures.

Here we describe an expanded set of cryptic lineages from multiple locations around the United States. While each sewershed contains its own signature lineages and at least some of the lineages appear to have diverged independently from one another, we present evidence that some likely shared a common ancestor. Finally, we show evidence of strong positive selection and rapid divergence of these lineages from ancestral SARS-CoV-2.

Results

Beginning in early 2021, wastewater surveillance programs including RBD amplicon sequencing (Fig 1A) were independently implemented in Missouri [9] and NYC [25]. A similar strategy was subsequently adopted in California. While the vast majority of sequences observed with this method matched to known lineages identified in patients, reproducible lineages that did not match the known circulating lineages were also detected. Herein, we refer to each RBD haplotype with a unique combination of amino acid changes as a lineage, and combinations of lineages that all have specific amino acid changes in common as lineage classes. Amino acid combinations identified that have not been seen previously from patients are referred to as cryptic lineages. Here we describe cryptic lineages detected from January 1, 2021 through March 15, 2022.

Lineage persistence and evolution over time

In total, cryptic lineages were observed in 9 sewersheds across 3 states (Table 1). Each cryptic lineage class was generally unique to a sewershed. These lineages contained between 4–24 non-synonymous substitutions, insertions, and deletions. In some cases, lineages were detected for a short duration but with multiple similar co-occurring sequences. For example, in Missouri sewershed MO33, a lineage class containing 4–5 RBD amino acid changes were consistently detected at low relative abundances from March 15 to the end of April 2021 (Fig. 1B, Table 1). A total of 7 unique sequences were spread across the 5 sampling events in this date range, and up to 5 unique sequences co-occurred within a given sample.

Table 1

Location	Date range when lineages appeared	Days within range	Number of samples	Number of RBD mutations
NY2	8/16/21–02/28/22	170	10	4–18
NY3	1/31/21 [8] −3/14/22	437	7	16–24
NY10	4/4/21–11/29/21	239	22	4–11
NY11	4/19/21–11/22/22	217	20	4–9
NY13	10/26/21–2/14/22	111	5	12–15
NY14	5/10/21–10/18/21	161	9	8–15
MO33	3/15/21–4/27/21	43	12	4–6
MO45	6/8/21–2/22/22	259	3	4–5
CA	11/4/21–12/21/21	47	3	16

Meanwhile, in other sewersheds, cryptic lineages were detected briefly, before disappearing, and then reappearing many months later. For example, in Missouri sewershed MO45, lineages were first detected in June 2021 and then were not seen again until February 2022 (Fig. 1C, Table 1). The longest observed lineage class was in sewershed NY3 where we previously reported a lineage class from January 2021 [8] that was detected sporadically until March 2022 (Fig. 2A, Table 1). On average, cryptic lineages lasted for around 6 months, such as the lineage class from NY14 which lasted from May to October, 2021 (Fig. 2B).

Fig. 2

NY3 and NY14 RBD amplifications. Overview of the SARS-COV-2 Spike RBD lineages identified from the A. NY3 and B. NY14 sewershed. Amino acid changes similar to or identical to changes in Omicron (BA.1) are indicated. The lineage previously referred to as WNY4 in [8] is indicated.

Each sewershed had its own unique set of lineages, but these lineages were not static. For instance, in NY10, the lineages first detected in April 2021 contained 4–5 RBD amino acid changes, but by October and November the lineages contained 8–11 RBD amino acid changes (Fig. 3, Table 1).

Fig. 3

NY10 RBD amplifications. Overview of the SARS-COV-2 Spike RBD lineages identified from the NY10 sewershed. Amino acid changes similar to or identical to changes in Omicron (BA.1) are indicated. Two lineages previously referred to as WNY1 and WNY2 in [8] are indicated.

In some cases, the sewersheds contained more than one lineage class. For instance, the NY11 sewershed contained several closely related lineages (class A) starting in April 2021, but a new set of lineages (class B) were detected starting in August 2021. These two classes were clearly distinct with very few amino acid changes in common (Fig. 4)

Fig. 4

NY11 RBD amplifications. Overview of the SARS-COV-2 Spike RBD lineages identified from the NY11 sewershed. Lineages designated A and B belong to two lineages groups that appear unrelated. Amino acid changes similar to or identical to changes in Omicron (BA.1) are indicated. The lineage previously referred to as WNY3 in [8] is indicated.

In addition to amino acid changes, several of the lineages observed in these sewersheds contained amino acid deletions near positions 445 and 484. For instance, lineages from NY2 contained a 444–445 deletion, NY11 and NY14 contained 443–444 deletion, NY3 and NY11 contained a deletion at position 484, and NY2 contained a deletion at position 483 (Fig. 2B, 4, 5A).

Fig. 5

NY2 and NY13 RBD amplifications. Overview of the SARS-COV-2 Spike RBD lineages identified from the A. NY3 and B. NY13 sewershed. Lineages designated A and B belong to two lineages groups that appear unrelated. Amino acid changes similar to or identical to changes in Omicron (BA.1) are indicated.

Most cryptic lineages detected did not appear to be derived from any known VOCs. The one exception was a lineage class containing amino acid changes N501Y and A570D in NY13 that first appeared on September 26, 2021, which suggested derivation from the Alpha VOC (Fig. 5B; Table 1). The Alpha VOC had been the dominant lineage in NYC between April and June 2021, but by September 26, 2021, it had been supplanted by Delta VOC and was no longer being detected in NYC [26].

Overall, specific lineage classes persisted within, but did not spread beyond, their individual sewersheds, with one notable exception. A cryptic lineage detected on August 16, 2021 in NYC sewershed NY2 precisely matched a lineage detected in sewershed NY11 between June-September 2021 (Fig. 4, 5A). While this precise lineage was never seen in NY2 again, several lineages with similar constellations of amino acid changes appeared in NY2 after October 4, 2021. The NY11 and NY2 sewersheds do not border each other, but are not separated by any bodies of water.

Rare and concerning amino acid changes are common in cryptic lineages and are sometimes shared with Omicron

In November 2021, the Omicron VOC was first detected in South Africa. This VOC contained eleven changes in the Spike protein between amino acids 410–510. Of these eleven amino acid changes, four (K417T, S477N, T478K, and N501Y) were present in previous VOCs. The remaining seven amino acid changes were rare prior to the Omicron VOC. All seven of these new amino acid changes had been detected in at least one of the wastewater lineages: N440K (MO33, NY11), G446S (NY2), E484A (MO45, NY10, NY11, NY2, NY13, CA), Q493R (NY3, NY14), G496S (NY2), Q498R (NY14, NY13), and Y505H (NY11, NY3, NY14, NY2, NY13, CA) (Figs. 1–6). None of the wastewater lineages have combinations of amino acid changes consistent with having a common ancestor with Omicron and most were initially detected prior to the emergence of Omicron. However, these shared amino acid changes suggest that the cryptic lineages were under selective pressures similar to those that shaped the Omicron lineage.

Fig. 6

Overview of the SARS-COV-2 Spike RBD lineages identified from the California sewershed. Amino acid changes similar to or identical to changes in Omicron (BA.1) are indicated.

Although each sewershed with cryptic lineages had its own signature combinations of amino acid changes, many of these changes were recurring among multiple sewersheds. Some of the more striking examples are described below.

N460K.

All nine of the sewersheds contained lineages with this change. Changes at this position are known to lead to evasion of class I neutralizing antibodies [27, 28]. However, this amino acid change is very rare, appearing in less than 0.01% of sequences in GISAID [22-24] submitted by March 15, 2022 (Table S1).

K417T.

Eight of the nine sewersheds contained lineages with the amino acid change K417T. Changes at this position are common and are known to participate in evasion from class I neutralizing antibodies [27, 28]. Although K417T was present in the Gamma VOC, K417N is the more common amino acid change at this position. The K417N amino acid change was not observed in any of the wastewater cryptic lineages.

N501S/T.

The amino acid changes N501S and N501T were seen in four and seven of the nine sewersheds, respectively. Changes at this position directly affect receptor binding and can affect the binding of multiple classes of neutralizing antibodies [19, 29, 30]. Although mutations at this position are very common, the most common change by far is N501Y, which was present in multiple VOCs. By contrast, N501S and N501T were present in less than 0.01% and 0.1% of sequences in GISAID submitted by March 15, 2022 (Table S1).

Q498H/Y.

Six of the nine sewersheds in this study contained lineages with the amino acid change Q498H or Q498Y. It should be noted that Q498Y differs from the Wuhan ancestral sequence by two nucleotide substitutions at the 498th codon (CAA→TAC). Q498H (CAA→CAC) is a necessary intermediary in this transition as TAA encodes a stop codon. In several cases both Q498H and Q498Y were seen in association with particular lineage classes including in NY2, 11, and 14 (Fig. 2B, 4, 5A) as well as a lineage class from California detected by the University of California, Berkeley wastewater monitoring laboratory (COVID-WEB) (Fig. 6). Changes at this position directly affect receptor binding and can affect the binding of multiple classes of neutralizing antibodies [19, 29, 30]. Notably, Q498H and Q498Y have been associated with mouse adapted SARS-CoV-2 lineages [31-33]. Both of these amino acid changes are very rare, appearing in less than 0.01% of sequences in GISAID submitted by March 15, 2022. Prior to November 2021, Q498Y had never been seen in a patient sample (Table S1).

E484A.

Six of the nine sewersheds contained lineages with the amino acid change E484A. Changes at this position are known to participate in evasion from class II neutralizing antibodies [27, 28]. Prior to the emergence of Omicron in November 2021, E484A was present in about 0.01% of sequences submitted to GISAID (Table S1).

Q493K.

Five of the nine sewersheds contained lineages with the amino acid change Q493K. Changes at this position directly affect receptor binding and can affect the binding of multiple classes of neutralizing antibodies [19, 27–30, 34]. This amino acid change is biophysically very similar to the Q493R mutation in Omicron. However, the Q493K amino acid change is very rare in patient derived sequences, appearing in less than 0.01% of sequences in GISAID submitted by March 15, 2022 (Table S1).

Y505H.

Five of the nine sewersheds contained lineages with the amino acid change Y505H. Prior to the emergence of Omicron in November 2021, Y505H was present in about 0.01% of sequences submitted to GISAID (Table S1).

K444T and K445A.

The amino acid changes K444T and K445A were each seen in four of the nine sewersheds. Changes at these positions are known to participate in evasion from class III neutralizing antibodies [28]. However, these amino acid changes are very rare, each appearing in less than 0.01% of sequences in GISAID submitted by March 15, 2022 (Table S1).

Y449R.

Three of the nine sewersheds contained lineages with the amino acid change Y449R. This change is noteworthy because, as of March 15, 2022, no sequences with this amino acid change had been submitted to GISAID (Table S1).

Long-read sequencing of S1 identifies substantial NTD modifications and suggests high dN/dS ratio

With each sample that contained novel cryptic lineages, attempts were made to amplify a larger fragment of the S1 domain of Spike. Amplification of larger fragments from wastewater is often inefficient, but sometimes can be achieved. To gain more information about the S1 domain of Spike and independently confirm the authenticity of the RBD lineages, we optimized a PCR strategy that amplifies 1.6 kb of the SARS-COV-2 Spike encompassing amino acids 57–579. These fragments were then either subcloned and sequenced or directly sequenced using Pacific Biosciences HiFi sequencing (Fig. 7A).

Fig. 7

S1 amplifications. A. Overview of the SARS-COV-2 Spike S1 lineages in the Alpha, Delta, Omicron VOCs and six of the sewersheds with cryptic lineages. S1 amplifications were sequenced by subcloning (SC) and Sanger sequencing, or were sequenced using a PacBio (PB) deep sequencing. B. Plot of the number of synonymous and non-synonymous changes in the S1 sequences shown.

The S1 amplification from the MO33 and MO45 sewersheds contained the RBD amino acid changes previously seen and each contained 3 additional amino acid changes upstream from the region sequenced using the targeted amplicon strategy described above (Fig. 7A). Many of the S1 amplifications from the NY10, NY11, NY13 and NY14 sewersheds contained numerous changes in S1 (Fig. 7A). In particular, many of the sequences contained deletions near amino acid positions 63–75, 144, and 245–248. All three of these areas are unstructured regions of the SARS-COV-2 spike where deletions have been commonly observed in sequences obtained from patients [35]. Two distinct S1 sequences were detected from the NY14 sample collected on June 28, 2021. Interestingly, the first sequence contained 13 amino acid changes which matched the RBD sequences from the same sewershed. The second sequence did not match any lineage that had been seen before, though it contained several mutations that were commonly seen in other cryptic lineages (see section 2.2). This second sequence presumably represented a unique lineage that had not been detected by routine surveillance.

A single S1 sequence was obtained from the NY13 samples collected on October 31, 2021. This sequence generally matched the RBD sequence from the same date, but did contain minor variations. Importantly, the S1 sequence contained deletions at positions 69–70 and 144, which, along with the amino acid changes N501Y and A570D, match the changes found in the Alpha VOC lineage. This information is consistent with the NY13 lineages being derived from the Alpha VOC.

Comparing the number of non-synonymous to synonymous mutations in a sequence can elucidate the strength of positive selection imposed on a sequence. The ratios of non-synonymous and synonymous mutations in this region of S1 from the Alpha, Delta, and Omicron VOCs (BA.1) were 19/0, 2/0, and 4/1, respectively. It was not possible to calculate the formal dN/dS ratios since many of the sequences did not have synonymous mutations in this region, so instead the numbers of non-synonymous and synonymous mutations were plotted. The cryptic lineages contained 5 to 25 total non-synonymous mutations and 0 to 2 total synonymous mutations (Fig. 7B).

Cryptic lineages from NCBI suggest an early common ancestor for many of the NYC lineages

In addition to RBD amplicon sequencing performed in our laboratories, we downloaded the 5609 SARS-CoV-2 wastewater fastq files from NCBI’s Sequence Read Archive (SRA) that were publicly available on NCBI on January 21, 2022 (not including submissions from our own groups). We screened these sequences for cryptic lineages by searching for recurring amino acid changes seen via RBD amplicon sequencing (K444T, Y449R, N460K, E484A, F486P, Q493K, Q493R, Q498H, Q498Y, N501S, N501T, and Y505H) (see above and Table S1), requiring at least two of these mutations with a depth of at least 4 reads. This strategy identified samples from 15 sewersheds (Table 2). Four were collected from unknown sewersheds in New Jersey and California in January 2021. The other 11 were collected by the company Biobot from NYC between June and August 2021. All but one of the lineages closely matched the cryptic sequences that had been observed via RBD amplicon sequencing from the same sewershed. The one exception was SRR16038150, which contained 4 amino acid changes that had not been seen in any of the previous sewershed samples in the same combination. The Biobot sequences were 40–96% complete and appeared to contain 30–100% cryptic lineages based on the frequency of mutation A23056C (Q498H/Y), a mutation shared with the lineages in all 11 sewershed samples from NYC. We speculate that the relative abundance of cryptic lineages was high because, during this period, NYC experienced the lowest levels of COVID-19 infections seen since the start of the pandemic. As a result, the sequences that matched the known circulating lineage were at low abundance.

Table 2.

Cryptic lineage whole genome sequences from nationwide surveys.

SRA Accession	State	Submitter	Sample Date	Percent cryptic lineage	Genome coverage	Sewershed	PANGO assignment	RBD Changes
SRR17120725	CA	Aquavitas	2021-01-04	7%	27,403	n/a[a]	NDb	E484A/Q498H/H519N
SRR16638981	NJ	Aquavitas	2021-01-18	7%	28,185	n/a[a]	NDb	E484A/Q498H/H519N
SRR16542155	NJ	Aquavitas	2021-01-18	7%	27,295	n/a[a]	NDb	E484A/Q498H/H519N
SRR16362183	NJ	Aquavitas	2021-01-04	100%	15,217	n/a[a]	NDb	E484A/Q498H/H519N
SRR16038150	NY	Biobot Analytics	2021-08-17	79%	28,227	NY2	B.1.503	Y449P/E484A/F490Y/Q498H
SRR16038156	NY	Biobot Analytics	2021-08-09	92%	24,595	NY11	B.1.503	K417T/K444T/Y449H/N460K/E484A/F490Y/Q498H
SRR15706711	NY	Biobot Analytics	2021-08-09	100%	11,877	NY11	NDb	K417T/K444T/Y449H/N460K/E484A/F490Y/Q498H/A570D
SRR15384049	NY	Biobot Analytics	2021-07-12	99%	24,001	NY10	B.1	Q493K/Q498Y/H519N/T572N)
SRR15291305	NY	Biobot Analytics	2021-07-05	100%	22,316	NY11	P.1.15	K417T/K444T/Y449H/E484A/F490Y/Q498H
SRR15291304	NY	Biobot Analytics	2021-07-04	100%	28,634	NY10	B.1	Q493K/Q498/H519N/T572N
SRR15202285	NY	Biobot Analytics	2021-06-28	100%	12,209	NY2	NDb	K444S/V445K/G446V/Y449R/L452Q/N460K/K462R/S477N/T478E/T478R/DEL483/E484P/F486I/F490P/G496S/Q498Y/P499S/N501T/Y505H/V511I
SRR15202284	NY	Biobot Analytics	2021-06-28	98%	16,281	NY14	NDb	K417T/K444S/DEL445-6/L452R/N460K/S477D/F486V/Q493K/Q498Y/P499S/N501T
SRR15202279	NY	Biobot Analytics	2021-06-28	30%	21,974	NY11	B.1	N440K/K444S/DEL445-6/L452Q/Y453F/N460K/S477N/D484/F486A/Q493K/Q498K/P499S/N501Y/H519N
SRR15128983	NY	Biobot Analytics	2021-06-16	99%	21,152	NY11	A.29	K444T/Y449H/E484A/Y489Y/F490Y/Q498H
SRR15128978	NY	Biobot Analytics	2021-06-16	100%	15,593	NY10	NDb	E484A/F486P/S494/Q498Y/H519N

n/a = not available;

ND = none designated

To compare the mutational profile among these different NYC samples, we first determined all of the mutations that occurred in at least 3 of the 11 cryptic lineages. We then produced a heat map to compare the frequency of each of these mutations from wastewater samples with the mutations that were reported from New York patient samples in June 2021 (Fig. 8). Surprisingly, the sewershed sequences often lacked two of the four consensus sequences that define the B.1 PANGO lineage (GISAID G clades or Nextstrain ‘20’ clades) of SARS-COV-2 [36]. Almost all patient samples collected in NYC during June 2020 contained the mutations C241T, C3037T, C14408T, and A23403G. The cryptic lineages from NYC wastewater all appeared to contain the mutations C3037T and A23403G, but possessed the ancestral sequences at positions 241 and 14408. In addition, there were two mutations in the S gene that were found in nearly all of the cryptic lineages, A23056C (Q498H/Y) and C24044T (L828F). Both of these mutations were found in less than 1% of patient samples. There were 3 additional mutations outside of the S gene that were highly prevalent in most of the wastewater samples, but essentially absent from patient samples: C25936G (Orf3 H182D), G25947C (Orf3 Q185H), and T27322C (Orf6 S41P). While other mutations were detected repeatedly within a sewershed, no other mutations spanned multiple sewersheds.

Fig. 8

Polymorphisms from wastewater genomes. Shown are all mutations present in at least three of the whole genome sequences from NYC listed in Table 2 and their corresponding amino acid changes. First column lists the prevalence of each mutation among all patients samples collected in June 2021 from New York. Each other column lists the prevalence of each mutation in each of the genome sequences.

To confirm that some of the cryptic lineages lacked the B.1 lineage consensus mutations, we designed primers to amplify and sequence the C14408 region of SARS-CoV-2 RNA isolated from wastewater. Indeed, samples from NY11 and NY10 that had a high prevalence of cryptic lineages were found to contain sequences that lacked C14408T (Fig. S1). However, when samples were amplified from the NY13 sewershed when the cryptic lineages there were present, we observed only the modern C14408T, as would be expected if the NY13 lineage were derived from the Alpha VOC. In addition, we performed whole genome sequencing on a March 30, 2021 sample from MO33 when the cryptic lineages were highly prevalent and did not detect any sequence that lacked C241T or C14408T, suggesting the cryptic lineages in this sewershed diverged after the emergence of the B.1 lineage (Fig. S2). Finally, we also analyzed the sequences from NCBI that contained the cryptic lineages from NJ and CA and did not find any sequences lacking C241T or C14408T. Thus, the lineages lacking C241T and C14408T appear to be limited to a subset of the cryptic lineages from NYC. It would appear that a SARS-CoV-2 lineage bearing mutations C3037T and A23403G, but possessing the ancestral genotype at positions 241 and 14408, was the direct ancestor of most of the cryptic lineages found in NYC.

Discussion

Our results point to the evolution of numerous SARS-CoV-2 lineages under positive immune selection whose source/host remains unknown.

Relatedness of and origin of cryptic lineages

We previously detected cryptic lineages via targeted amplicon sequencing [8], but lacked information about their derivation. Here, from comparison of the sewersheds for which whole genome sequencing is available, it is clear that the cryptic lineages from wastewater are not all derived from a common ancestor. The NY13 lineage appeared to be derived from the Alpha VOC. If this is true, the NY13 lineage most likely branched off from Alpha sometime in early to mid-2021 when that variant was common in NYC. However, many lineages from the NY10, NY11, NY2, and NY14 sewersheds in New York appear to share a common ancestor that branched off from a pre-B.1 lineage. Additionally, we often observed swarms of related sequences that co-occurred within a sewershed on a single date, and accumulated new mutations over time, suggesting continued diversification from a single origin within each sewershed.

Comparison with the Omicron VOC

The Omicron VOC and the wastewater lineages appear to have been subjected to high positive selection. While prior VOCs had 3 or fewer amino acid changes in the amplified region of the RBD, the Omicron VOC (BA.1) contained 11 and the cryptic lineages from wastewater averaged over 10. By comparison, a cluster of SARS-COV-2 sequences that appear to have circulated in white-tailed deer for over a year accumulated only 2 amino acid changes in this region [37]. Of the nonsynonymous RBD mutations in Omicron, four were in at least one prior VOC: K417N, S477N, T478K, and N501Y. The other seven were relatively rare; N440K was present in 0.2% of sequences and the other six were each present in less than 0.1% of sequences in GISAID prior to November 1, 2021. All of the rare Omicron changes were observed in at least one of the cryptic wastewater lineages. Collectively, this suggests that the wastewater lineages and the Omicron VOC likely arose under similar selective pressures. The high dN/dS ratios found in cryptic lineages and in Omicron suggest that these selective pressures must be exceptionally strong.

Source of Lineages

In spite of detailed tracking and cataloging of the cryptic lineages, the question where they are coming from remains unanswered. The most parsimonious explanations are 1) undetected spread within the human population, 2) prolonged shedding by individuals, most likely immunosuppressed, or 3) spread in animal reservoirs.

Undetected spread in the population appears unlikely. While the sequencing rate for US patient samples is not 100%, it is high enough that population-level spread of cryptic lineages would not be missed. Alternatively, as it is known that SARS-CoV-2 can replicate in gastrointestinal sites [38, 39], the lack of detection of cryptic lineages by clinical sequencing could be explained by the potential adaptation of some SARS-CoV-2 to replicate exclusively in the gastrointestinal tract [1, 38]. Nonetheless, even if replication of these lineages were occurring outside of the nasopharyngeal region, this could not explain why cryptic lineages generally remain geographically constrained.

The simplest explanation for the appearance of cryptic lineages in wastewater is that they are shed by immunosuppressed patients with persistent infections. Indeed, the vast majority of amino acid changes in the RBD of the Omicron VOC and the cryptic lineages confer resistance to neutralizing antibodies. In particular, substitutions at positions 417, 440, 460, 484, 493, 498 and 501 have all been well documented to lead to immune evasion [17, 27, 34, 40–42]. Additionally, RBD changes K417T, N440K, N460K, E484A, Q493K, and N501Y have all been observed in persistent infections of immunocompromised patients [43, 44]. Given the repeated appearance of these mutations in diverse sewersheds, the majority of the selective pressure on the cryptic lineages is almost certainly immune pressure. The counterargument to this explanation is the sheer volume of viral shedding required to account for the wastewater signal. Many of the sewersheds process 50–100 million gallons of wastewater per day. Reliable amplification of a sequence from wastewater generally requires that the sequence is present at least 10,000 copies per liter. Therefore, detection of a specific virus lineage in such a sewershed would seem to require several trillion virus particles to be deposited each day. If this signal were derived from a single infected patient or even a small group of patients, those patients would have to shed exponentially more viruses than typical COVID-19 patients.

The final explanation for the cryptic lineages in wastewater is that they are shed into wastewater by an animal host population. Previously, we determined through rRNA analysis of several NYC sewersheds that the major non-human mammals that contribute to the wastewater are cats, rats, and dogs [8]. Of these three, rats were the only species that seemed to be a plausible candidate. Indeed, we also showed that the cryptic lineages from the sewersheds had the ability to utilize rat and mouse ACE2 [8]. However, one of the sewersheds with the most consistent signal in 2021 was NY10, which had little to no rat rRNA. In addition, it is not clear why circulation in an immune competent animal such as a dog or a rat would result in a more rapid selection of immune escape mutations than circulation in humans, yet the cryptic lineages display accumulation of many times more immune escape changes than seen in viruses circulating in the human population.

The importance of wastewater sequencing methodology for identification of novel variants

In order to provide information regarding the appearance and spread of SARS-CoV-2 variants in communities, next generation sequencing technologies have been applied to sequence SARS-CoV-2 genetic material obtained from sewersheds around the world [45-47]. Commonly, SARS-CoV-2 RNA extracted from wastewater is amplified using SARS-CoV-2 specific primers that cover the entire genome [48-50]. Bioinformatic pipelines are employed to identify circulating SARS-CoV-2 variants [16, 51]. In general, the presence and abundance of variants in wastewater corresponds to data obtained from clinical sequencing [45, 46]. However, to our knowledge, there have been no other reports of cryptic lineages detected in wastewater that were not also observed in clinical sequence data. A major issue with generating whole genome sequence data from nucleic acid isolated from wastewater is sequence dropout over diagnostically important regions of the genome [48, 52, 53]. In some cases, diagnostically important regions of the genome that accumulate many mutations, such as the Spike RBD, receive little to no sequence coverage, making variant attribution difficult. Since wastewater contains a mixture of virus lineages and whole genome sequencing relies on sequencing of small genome fragments, mutations appearing on different reads cannot be linked together. Indeed, some variant identification pipelines map reads to reference genomes to estimate the probability that mutations are found in the same genome [16]. Such strategies would not be able to detect variants containing unique constellations of mutations. Detecting novel variants that are present at low relative abundances may be better achieved by targeted amplicon sequencing, such as the strategy we present here.

Summary

Over the past 15 months, cryptic SARS-CoV-2 lineages never seen in human patients have appeared in community wastewater in several locations across the USA [8]. These lineages have persisted, intermittently, often as swarms of closely related haplotypes that acquired additional amino acid changes over time, for up to 14 months. Evidence suggests that some of the lineages may have arisen during the initial phases of the pandemic in early to mid-2020. Significantly, these lineages often contained amino acid changes that have rarely or never appeared in contemporaneous variants, at least until the appearance of the Omicron VOC. Many of these amino acid changes are associated with evasion of antibody-mediated neutralization. Collectively, nonsynonymous substitutions in these lineages overwhelmingly outnumbered synonymous substitutions, indicating that these lineages have undergone exceptionally strong positive selection.

Three hypotheses for the origins of these lineages have been proposed: 1) undetected transmission, 2) long-term infections of immunocompromised patients and 3) possible animal reservoirs. Although immunosuppressed populations are the simplest explanation, it is difficult to reconcile the magnitude of the signal with individual patients being the source. Regardless of the origins and dynamics of cryptic variant shedding, our results highlight the ability of wastewater-based epidemiology to more completely monitor SARS-CoV-2 transmission and genetic diversity than can patient based sampling, at scale and at a greatly reduced cost. Given that multiple VOCs may have gone undetected until suddenly appearing, highly mutated, in apparently single evolutionary leaps [12], it is crucial to the early detection of the next variant of concern that novel SARS-CoV-2 genotypes are monitored for evidence of significant expansion. Importantly, patient sampling efforts, despite occurring with an intensity not seen in any prior epidemic, were unable to identify intermediary forms of most VOCs. Monitoring of wastewater, particularly using a targeted sequencing approach, likely provides the best avenue for detecting developing VOCs.

Materials and Methods

Wastewater sample processing and RNA extraction

24-hr composite samples of wastewater were collected weekly from the inflow at each of the wastewater treatment plans.

NYC:

Samples were processed on the day they were collected and RNA was isolated according to our previously published protocol [6]. Briefly, 250 mL from a 24-hr composite wastewater sample from each WWTP were centrifuged at 5,000 × g for 10 min at 4°C to pellet solids. A 40 mL aliquot from the centrifuged samples was passed through a 0.22 μM filter (Millipore). To each corresponding filtrate, 0.9 g sodium chloride and 4.0 g PEG 8000 (Fisher Scientific) were added. The tubes were kept at 4°C for 24 hrs and then centrifuged at 12,000 × g for 120 minutes at 4 °C to pellet the precipitate. The pellet was resuspended in 1.5 mL TRIzol (Fisher Scientific), and RNA was purified according to the manufacturer’s instructions.

MO:

Samples were processed as previously described [9]. Briefly, wastewater samples were centrifuged at 3000×g for 10 min and then filtered through a 0.22 μM polyethersolfone membrane (Millipore, Burlington, MA, USA). Approximately 37.5 mL of wastewater was mixed with 12.5 mL solution containing 50% (w/vol) polyethylene glycol 8000 and 1.2 M NaCl, mixed, and incubated at 4°C for at least 1 h. Samples were then centrifuged at 12,000×g for 2 h at 4°C. Supernatant was decanted and RNA was extracted from the remaining pellet (usually not visible) with the QIAamp Viral RNA Mini Kit (Qiagen, Germantown, MD, USA) using the manufacturer’s instructions. RNA was extracted in a final volume of 60 μL.

CA:

Samples were processed as previously described [54]. Briefly, 40 mLs of influent was mixed with 9.35g NaCl and 400 uL of 1M Tris pH 7.2, 100mM EDTA. Solution was filtered through a 5-um PVDF filter and 40 mLs of 70% EtNY11 was added. Mixture was passed through a silica spin column. Columns were washed with 5 mL of wash buffer 1 (1.5 M NaCl, 10 mM Tris pH 7.2, 20% EtNY11), and then 10 mL of wash buffer 2 (100 mM NaCl, 10 mM Tris pH 7.2, 80% EtNY11). RNA was eluted with 200 ul of ZymoPURE elution buffer.

Targeted PCR: MiSeq sequencing

The primary RBD RT-PCR was performed using the Superscript IV One-Step RT-PCR System (Thermo Fisher Scientific, 12594100).  Primary RT-PCR amplification was performed as follows: 25 °C (2:00) + 50 °C (20:00) + 95 °C (2:00) + [95 °C (0:15) + 55 °C (0:30) + 72 °C (1:00)] × 25 cycles using the MiSeq primary PCR primers CTGCTTTACTAATGTCTATGCAGATTC and NCCTGATAAAGAACAGCAACCT. Secondary PCR (25 μL) was performed on RBD amplifications using 5 μL of the primary PCR as template with MiSeq nested gene specific primers containing 5′ adapter sequences (0.5 μM each) acactctttccctacacgacgctcttccgatctGTRATGAAGTCAGMCAAATYGC and gtgactggagttcagacgtgtgctcttccgatctATGTCAAGAATCTCAAGTGTCTG, dNTPs (100 μM each) (New England Biolabs, N0447L) and Q5 DNA polymerase (New England Biolabs, M0541S).  Secondary PCR amplification was performed as follows: 95 °C (2:00) + [95 °C (0:15) + 55 °C (0:30) + 72 °C (1:00)] × 20 cycles. A tertiary PCR (50 μL) was performed to add adapter sequences required for Illumina cluster generation with forward and reverse primers (0.2 μM each), dNTPs (200 μM each) (New England Biolabs, N0447L) and Phusion High-Fidelity or (KAPA HiFi for CA samples) DNA Polymerase (1U) (New England Biolabs, M0530L).  PCR amplification was performed as follows: 98 °C (3:00) + [98 °C (0:15) + 50 °C (0:30) + 72 °C (0:30)] × 7 cycles +72 °C (7:00).  Amplified product (10 μl) from each PCR reaction is combined and thoroughly mixed to make a single pool. Pooled amplicons were purified by addition of Axygen AxyPrep MagPCR Clean-up beads (Axygen, MAG-PCR-CL-50) or in a 1.0 ratio to purify final amplicons. The final amplicon library pool was evaluated using the Agilent Fragment Analyzer automated electrophoresis system, quantified using the Qubit HS dsDNA assay (Invitrogen), and diluted according to Illumina’s standard protocol. The Illumina MiSeq instrument was used to generate paired-end 300 base pair reads. Adapter sequences were trimmed from output sequences using Cutadapt.

Long PCR and subcloning.

The long RBD RT-PCR was performed using the Superscript IV One-Step RT-PCR System (Thermo Fisher Scientific, 12594100). Primary long RT-PCR amplification was performed as follows: 25 °C (2:00) + 50 °C (20:00) + 95 °C (2:00) + [95 °C (0:15) + 55 °C (0:30) + 72 °C (1:30)] × 25 cycles using primary primers CCCTGCATACACTAATTCTTTCAC and TCCTGATAAAGAACAGCAACCT. Secondary PCR (25 μL) was performed on RBD amplifications using 5 μL of the primary PCR as template with nested primers (0.5 μM each) CATTCAACTCAGGACTTGTTCTT and ATGTCAAGAATCTCAAGTGTCTG, dNTPs (100 μM each) (New England Biolabs, N0447L) and Q5 High-Fidelity DNA Polymerase (New England Biolabs, M0491L). Secondary PCR amplification was performed as follows: 95 °C (2:00) + [95 °C (0:15) + 55 °C (0:30) + 72 °C (1:30)] × 20 cycles.

Positive amplifications were visualized in an agarose gel stained with ethidium bromide, excised, and purified with a NuceloSpin Gel and PCR Clean-up Kit (Macherey-Nagel, 74609.250). Gel purified DNA was subcloned using a Zero Blunt TOPO PCR Cloning Kit (Invitrogen, K2800–20SC). Individual colonies were transferred to capped test tubes containing 10 ml of 2X YT broth (ThermoFisher, BP9743–5). Test tubes were incubated at 37°C and shook at 250 rpm for 24 hours. The resulting E. Coli colonies were centrifuged for 10 minutes at 5000 × g and the supernatant was decanted. Plasmid DNA was extracted from the pellet using a GeneJet Plasmid Miniprep Kit (ThermoFisher, K0503). The concentration of plasmid DNA extracts was measured using a NanoDrop One (ThermoFisher, ND-ONE-W).

Subcloning

The 1.6kb S1 fragment was cloned into the pMiniT vector using the NEB PCR Cloning Kit (NEB #E1202) protocol. Briefly, 5 ul of the RT-PCRed fragment was ligated to the linearized pMiniT 2.0 Vector and transformed into NEB 10-beta Competent E. coli (NEB #C3019). Transformed cells were outgrown in 950ul of NEB 10-beta/Stable Outgrowth Medium at 37°C and shaken at 200 rpm for 1 hour. Dilutions of 1:10, 1:100, and 1:1000 were plated onto 100 ug/mL LB-ampicillin using glass beads. Plates were incubated overnight at 37°C. Following this, single colonies were swatched onto 100ug/mL LB-ampicillin plates and pools of 5 colonies were lysed in 20ul H20 for 5 minutes at 95°C. 5 ul of lysed template was added to the following PCR mix: 10 ul Q5 High-Fidelity 2X Master Mix (NEB #M0492S) + 0.8 ul 5uM CATTCAACTCAGGACTTGTTCTT forward primer (2402F) + 0.8 ul 5uM ATGTCAAGAATCTCAAGTGTCTG (2376R) reverse primer+ 5 ul H2O. Thermocycling conditions were as follows: 95°C (2:00) + [95°C (00:15) + 55°C (00:30) + 72°C (3:00)] × 40 cycles + 72°C (1:00). PCR products were run on a 0.8% agarose gel and successfully cloned colonies were determined after ethidium bromide staining. Positive pools were then sent for Genewiz Sanger sequencing with the following primer pairs: 2402F (5uM), 2376R (5uM), and NEB Cloning Kit Cloning Analysis Forward Primer (100uM), Cloning Analysis Reverse Primer (100uM). The following additional primer pairs were included as needed: TGCGAATAATTGCACTTTTGA and TGCTACCGGCCTGATAGATT, GGACCTTGAAGGAAAACAGG and TGCTACCGGCCTGATAGATT, and ATCTCCCTCAGGGTTTTTCG and CCATTACAAGGTGTGCTACCG. Pools with mixed Sanger sequencing signals were resequenced as individual colonies picked from the swatch plates and DNA was isolated using the QIAprep Spin Miniprep Kit (#27106).

PacBio sequencing

A nested RT-PCR protocol was used to generate 1.6kb Spike amplicons from wastewater RNAs for PacBio sequencing. The primary RT-PCR amplification was performed with the SuperscriptTM IV One-Step RT-PCR System (Invitrogen) and the same thermal cycling program as described above for MiSeq amplicons. These inter Spike gene-specific primer sequences (5’-[BC10ab]-ATTCAACTCAGGACTTGTTCTT and 5’-[BC10xy]-ATGTCAAGAATCTCAAGTGTCTG) were tagged directly on their 5’ ends with standard 16 bp PacBio barcode sequences (Supplemental Table Z) and used with asymmetric barcode combinations that allow large numbers of samples to be pooled prior to sequencing. The following thermal cycling profile was used for nested PCR: 98 °C (2 min) + [98 °C (10 sec) + 55 °C(10 sec) + 72 °C (1 min)] × 20 cycles + 72 °C (5 min). The resulting PCR amplicons were then subjected to three rounds of purification with AMPure XP beads (Beckman Coulter Life Sciences) in a ratio of 0.7:1 beads to PCR. Purified amplicons were quantified using a Qubit™ dsDNA HS kit (ThermoFisher Scientific) and pooled prior to PacBio library preparation.

After ligation of SMRTbell adaptors according to the manufacturer’s protocol, sequencing was completed on a PacBio Sequel II instrument (PacBio, Menlo Park, CA USA) in the Genomic Sequencing Laboratory at the Centers of Disease Control in Atlanta, GA, USA. Raw sequence data was processed using the SMRT Link v10.2 command line toolset (Software downloads - PacBio). Circular consensus sequences were demultiplexed based on the asymmetric barcode combinations and subjected to PB Amplicon Analysis to obtain high-quality consensus sequences and search for minor sequence variants.

Bioinformatics

MiSeq and PacBio processing

Sequencing reads were processed as previously described. Briefly, VSEARCH tools were used to merge paired reads and dereplicate sequences [55]. Dereplicated sequences from RBD amplicons were mapped to the reference sequence of SARS-CoV-2 (NC_045512.2) spike ORF using Minimap2 [56]. Mapped amplicon sequences were then processed with SAM Refiner using the same spike sequence as a reference and the command line parameters “--Alpha 1.8 --foldab 0.6” [9].

The covariant deconvolution outputs were used to generate the haplotype plots in figures 1–7. MiSeq sequences from those outputs were collected by sewershed and multiple runs of the same sample averaged. The collected MiSeq sequences were processed to remove the Alpha, Beta, Gamma, Delta and Omicron VOC lineage sequences and remove individual polymorphisms that appeared only once in sewersheds that had unknown lineages on more than 3 dates. Sequences that appeared in only one sample and had an abundance less than .02 were discarded. The resulting sequences were then plotted along with their date and abundance. The PacBio sequences were similarly collected to generate the haplotype plot in Fig. 7. The covariant outputs, their collection and variant lineage processing outputs, and the scripts for processing them are available upon request. The MO raw sequence reads are available in NCBI’s SRA under the BioProject accession PRJNA748354. The NY raw sequence reads are available in NCBI’s SRA under the BioProject accession PRJNA715712.

NCBI SRA screening

Raw reads were downloaded and then processed similar to MiSeq sequencing except the reads were mapped to the entire SARS-CoV-2 genome and SAM Refiner was run with the parameters ‘--wgs 1 --collect 0 --indel 0 --covar 0 --min_count 1 --min_samp_abund 0 --min_col_abund 0 --ntabund 0 --ntcover 1’. Unique sequence outputs from SAM Refiner were then screened for specific amino acid changes with a custom script (available upon request). The nt call outputs of samples of interest were used to determine other variations in the genomes sequenced.

14408 sequencing

The long RBD RT-PCR was performed using the Superscript IV One-Step RT-PCR System (Thermo Fisher Scientific, 12594100). Primary long RT-PCR amplification was performed as follows: 25 °C (2:00) + 50 °C (20:00) + 95 °C (2:00) + [95 °C (0:15) + 55 °C (0:30) + 72 °C (1:30)] × 25 cycles using primary primers ATACAAACCACGCCAGGTAG and AACCCTTAGACACAGCAAAGT. Secondary PCR (25 μL) was performed on RBD amplifications using 5 μL of the primary PCR as template with nested primers (0.5 μM each) ACACTCTTTCCCTACACGACGCTCTTCCGATCTGGTAGTGGAGTTCCTGTTGTAG and GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTAGCACGTAGTGCGTTTATCT, dNTPs (100 μM each) (New England Biolabs, N0447L) and Q5 High-Fidelity DNA Polymerase (New England Biolabs, M0491L). Secondary PCR amplification was performed as follows: 95 °C (2:00) + [95 °C (0:15) + 55 °C (0:30) + 72 °C (1:30)] × 20 cycles.

Whole genome sequencing

Whole genome sequencing of the SARS-CoV-2 genome from the MO33 sewershed was performed using the NEBNext ARTIC SARS-CoV-2 Library Prep Kit (Illumina). Amplicons were sequenced on an Illumina MiSeq instrument. Output sequences were analyzed using the program SAM Refiner [57].