Monitoring SARS-CoV-2 Circulation and Diversity through Community Wastewater Sequencing, the Netherlands and Belgium

1Since its discovery, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has caused >100 million confirmed cases of coronavirus disease (COVID-19). The global effects of SARS-CoV-2 and the need to learn more about its origin and epidemiology have resulted in the sequencing of >416,000 genomes as of January 2021. This work has enabled the identification of groups of viruses that, on the basis of their genetic diversity, can be associated with geographic and temporal patterns of virus spread. Nextstrain (https://nextstrain.org) currently divides SARS-CoV-2 diversity into 12 major global clades (19A, 19B, and 20A–20J), on the basis of high prevalence, signature mutations, and geographic spread.

2Although SARS-CoV-2 primarily affects respiratory tract tissues, it can also replicate in the gastrointestinal tract, as evidenced by in vitro infection of enteroids, presence of viral proteins in gastrointestinal epithelium biopsy specimens, and detection of infectious virus in stool samples. Viral RNA is shed in the feces of ≈40% of infected persons, often for longer periods than the virus can be detected in nasal swab specimens. SARS-CoV-2 RNA has been detected in urine occasionally (<5% of infected patients).

3Because of the rapid spread of SARS-CoV-2, individual screening of clinical cases and study of viral diversity on a population level are challenging. Various reports have demonstrated that enteric and respiratory viruses can be detected in wastewater. This finding has led to the recognition of wastewater-based epidemiology as a potentially valuable tool to assess the spread of the disease at a community level. Recently, the Water Research Institute in the Netherlands and other groups have demonstrated temporal correlations between SARS-CoV-2 RNA titers in sewage and the number of reported cases in a city or county when >26 gene copies per liter could be detected. Therefore, sewage testing is currently considered globally to be an adjunct to patient-based surveillance and demonstrates promise as an early warning indicator of increasing virus circulation.

4Enhanced surveillance is a key pillar of the current strategy to control the spread of SARS-CoV-2 and includes frequently testing mildly symptomatic persons, investigating infection clusters to identify possible common exposures, and monitoring hospital admission rates. Whole-genome sequencing of SARS-CoV-2 from clinical samples has been adopted as an additional tool to identify clusters. Particularly in geographic areas with minimal virus circulation, sequencing can help identify possible sources, provided that sufficient background sequencing has been performed. So far, little work has been done to correlate SARS-CoV-2 diversity in sewage samples with diversity in patients. We used next-generation sequencing (NGS) of SARS-CoV-2 from wastewater samples to assess whether these samples reflect the diversity of SARS-CoV-2 circulating within the population of the Netherlands and Belgium.

5Wastewater specimens were collected as 24-h flow-dependent composite samples and processed as previously described. Debris of 100–200 mL of sewage samples was pelleted and the supernatant was concentrated by using 100 kDa Centricon ultrafilters (Millipore Sigma, https://www.emdmillipore.com); in vitro–transcribed dengue virus type-2 RNA was added as an internal extraction control. RNA was extracted by using the Nuclisens kit (bioMérieux, https://www.biomerieux.com) and KingFisher purification system (Thermo Fisher Scientific, https://www.thermofisher.com). RNA was screened by quantitative reverse transcription PCR (qRT-PCR) with 5 primer–probe sets targeting the SARS-CoV-2 nucleocapsid (N) gene (N1–N3), envelope (E) gene for all sarbecoviruses, and the internal control.

6We performed SARS-CoV-2–specific multiplex PCR for nanopore sequencing as described previously. Primers for 89 overlapping amplicons spanning the genome were used in 2 PCR pools. Libraries were generated by using the Oxford Nanopore native barcode kits (Oxford Nanopore Technologies, https://nanoporetech.com) and sequenced on a R9.4 flow cell.

7Illumina sequencing was performed as described previously. Amplicons were generated by the multiplex PCR described previously. Amplicons were purified with 0.8X AMPure XP beads (Beckman Coulter, https://www.beckmancoulter.com) and 100 ng of DNA was converted into paired-end Illumina sequencing libraries by using the KAPA HyperPlus library preparation kit (Roche, https://www.roche.com). We used the KAPA Unique Dual-Indexed Adapters Kit (Roche) to enable subsequent sequencing of multiple libraries in a single Illumina MiSeq version 3 flowcell (2 × 300 cycles) (Illumina, https://www.illumina.com).

8Raw sequence data were processed as previously described. We used a snakemake script to demultiplex fastq raw reads by using Porechop (https://github.com/rrwick/Porechop), to trim primers by using Cutadapt, and to perform a reference-based alignment by using minimap2 to GISAID sequence EPI_ISL_412973 (https://www.gisaid.org). The run was monitored by using RAMPART (https://artic-network.github.io/rampart). The consensus genome was extracted by using 2 analyses for which positions with a coverage <10X or <30X were replaced with an N. We confirmed mutations in the genome by manually checking the alignment in Ugene and resolved homopolymeric regions by consulting reference genomes. On the basis of previous studies, we considered mutations with >30X coverage high quality, whereas mutations >10X and <30X coverage were considered low quality.

9We used a customized Galaxy workflow for all processing, reference-based alignment, and variant analysis. Raw sequencing reads were filtered by using Fastp to remove adaptor contamination, ambiguous bases, low quality reads (Phred score <30), and fragments <50 nt. Reads were mapped against GISAID sequence EPI_ISL_412973 by using the default settings of BWA-MEM (H. Li, unpub. data, https://arxiv.org/abs/1303.3997). Reads were realigned by using the leftalign utility from FreeBayes (E. Garrison, unpub. data, https://arxiv.org/abs/1207.3907). All reads with mapping scores of <30 were discarded. Consensus sequences and variants were generated by using iVar. Final consensus sequences (frequency >50%) were constructed by using all mapped reads with a coverage of >5X and Phred score of >30. For detection of low-frequency variants (LFVs), we used parameters as follows: minimum coverage of 50X, Phred score >30, and a minimum frequency threshold of 10%. Variant calling was confirmed by manual inspection of the aligned reads in Ugene. Variant positions are given with respect to the Wuhan-Hu-1 strain (MN908947). We uploaded all consensus sequences with coverage >50% to GISAID (accession nos. EPI_ISL_539300–25).

10The first dataset included all full-length SARS-CoV-2 genomes from the Netherlands (1,544 genomes) and Belgium (888 genomes) from GISAID as of July 8, 2020. The second dataset was a subsample representative of the global diversity of all SARS-CoV-2 sequences in GISAID as of June 1, 2020. This global dataset contained 2,552 subsampled sequences (full length with Ns <5%) to include 1 unique genome per country or state per week. We aligned sequences with >75% genome coverage by using MAFFT (https://mafft.cbrc.jp/alignment/server) and inferred maximum-likelihood trees by using the best predicted models general time-reversible plus F plus R3 (global subsample) and general time-reversible plus F plus R2 (Netherlands–Belgium dataset) and bootstrap with 1,000 replicates. Trees were visualized by using Figtree version 1.4.4 (http://tree.bio.ed.ac.uk/software/figtree). Clades were assigned by using the Nextclade tool.

11Previously, sewage samples collected from 6 locations in the Netherlands (and Schiphol Airport) were tested by qRT-PCR to investigate the levels of SARS-CoV-2 RNA. To further investigate the genetic diversity of SARS-CoV-2, we subjected 55 wastewater samples obtained from 13 locations in the Netherlands (48 samples) and 7 locations in Belgium (7 samples) with cycle threshold (Ct) values of <36 to whole-genome sequencing by using nanopore technology. The wastewater treatment plants in the Netherlands served ≈200,000–980,000 inhabitants; Schiphol was estimated to serve 54,000 persons. The samples covered a period of 70 days (March 25–June 3, 2020); of all 55 samples, 2 (Franeker-92719 and AmsterdamWest-92852) were sequenced by nanopore twice. Of the 55 samples, 24 were also sequenced by Illumina (Table 1).

12We used 4 primer–probe sets targeting the N (N1–N3) genes and E gene to evaluate the concentration of SARS-CoV-2 in sewage samples (Table 1). The percentage of the genome covered by the assembly of nanopore reads (>10X coverage) ranged from 0% to 99.2%. We found an inverse sigmoidal correlation between the percentage of the genome assembled from nanopore sequencing reads and the N and E gene Ct values (Figure 1). The Ct values at which half of the genome could be obtained were 34.6 for N1, 33.8 for N2, 33.2 for N3, and 32.5 for E. No correlation was observed between Ct values and the percentage of the genome assembled from Illumina reads (Appendix Figure 1).

13We performed phylogenetic analysis to assess whether consensus sequences from sewage could be associated with clinical samples from the same region. A total of 22 genomes (20 from nanopore and 2 from Illumina runs) with a coverage >75% of the genome were obtained from 20 samples. We used these sequences to infer a maximum-likelihood tree using all sequences from the Netherlands and Belgium available in GISAID and a maximum-likelihood tree using a subset representative of the global diversity of SARS-CoV-2 in GISAID. In general, the sequences from the Netherlands and Belgium grouped into 5 clades (Figure 2, panel A), and most of the sequences belonged to clade 20A (52.0% for the Netherlands and 47.7% for Belgium). The clades 19B and 20C were less prevalent; 8.9% of sequences from the Netherlands belonged to 19B and 1.2% to 20C, whereas 10.4% of Belgium sequences belonged to 19B and 0.3% to 20C. Both trees showed that sewage samples grouped within clades 19A, 20A, and 20B (Figure 2). Samples Franeker-92719 and HeeswijkDinther-92499 clustered with sequences isolated from patients from the same region (Figure 2, panel A), indicating that sewage samples can be linked to specific outbreaks. Included in the phylogenetic trees were 2 samples with 2 consensus sequences (AmsterdamWest-92852 and Franeker-92719), which demonstrated 2-mutation differences between consensus sequences of the same sample (Appendix Table 1). Despite this discrepancy, consensus sequences from the same sample clustered within the same clade (Appendix Figures 2, 3). Some sequences clustered close to the root of the tree, probably because of the presence of multiple strains within 1 sample, which resulted in a combination of mutations in their consensus sequences.

14To associate samples with a particular clade or cluster, we compared all consensus sequences, including partial sequences, with the Wuhan-Hu-1 reference isolate. A total of 145 single-nucleotide polymorphisms (SNPs) were detected in our dataset (Appendix Table 1). Of these, 24 SNPs were detected in >1 sequence. We also detected SNPs in the Netherlands sewage sequences with a geographic regional signal, which were present in the Netherlands clinical samples at much higher frequencies than in global or Belgium clinical samples, such as T514C and C1594T (Appendix Table 2).

15Finding clade-defining mutations in the consensus sequence suggests the dominance of a certain clade within a sample; the presence of these mutations can also aid in the detection of virus mixtures in a sample. During the period of wastewater-sample collection, Nextstrain defined 5 major clades (19A, 19B, 20A, 20B, and 20C). Each clade is defined by the presence of >2 linked mutations. Clade 19A is the root clade and contains the Wuhan-Hu-1 reference sequence. Both 19B and 20A emerged from 19A, where 2 and 3 linked mutations define these major clades: T28144C and C8782T define 19B; and C3037T, C14408T, and A23403G define 20A. Clades 20B and 20C emerged from 20A, where the trinucleotide substitution GGG28881–28883AAC defines 20B and the linked mutations C1059T and G25563T define 20C. Nucleotide substitution A23403G, a signature mutation of clades 20A, 20B, and 20C that generates the D614G amino-acid substitution in the S glycoprotein, was detected in 83.6% (51/61) of the samples that were sequenced at this region (Appendix Table 1). The GGG28881–28883AAC substitution was detected in 41.9% (18/43) of the sequences. One of the 2 mutations defining the low-prevalence clades 20C and 19B (C1059T and T28144C) was found in 2 and 3 consensus sequences. However, these sequences could not be assigned to these clades because regions containing the additional clade-defining mutations were not sequenced with sufficient coverage. The hCoV-19/env/Netherlands/Amersfoort-92503-N/2020 sequence contained a mix of clade-defining mutations: C1059T, which defines 20C; T28144C, which defines 19B; and GGG28881–28883AAC, which defines 20B. This finding indicates that the obtained consensus sequence does not represent a single strain.

16In addition to the clade-defining mutations, we detected 49 and 63 SNPs that were not present in either the Netherlands (1,544 sequences) or Belgium (888 sequences) datasets but were seen in the global dataset (55,074 sequences), although with <1% prevalence (Appendix Table 2). Moreover, we detected 51 novel mutations in sewage consensus sequences that were not previously reported, of which 48 were supported by coverage above the thresholds set for high quality (coverage >30× for Nanopore and coverage >5× and Phred score >30 for Illumina). Discrepancies between consensus sequences of the same sewage sample can occur. AmsterdamWest-92852 was sequenced 3 times and 4 positions varied (Appendix Table 1). These differences are explained by the presence of variant sites in a single sample in similar percentages, which resulted in differences in consensus sequences between sequencing runs.

17Given that sewage samples are likely to contain a mixture of SARS-CoV-2 strains, we performed a variant analysis with Illumina data to distinguish multiple strains within single samples. By using a coverage >50×, Phred score >30, and a frequency threshold of >10% as settings, we found 21 positions with at least 1 sample containing major and minor variants (Table 2). Of these, 14 mutations resulted in changes at the amino acid level (12 nonsynonymous mutations and 2 deletions). Of note, 8 of these (4497C, 10514C, 11484T, 13046A, 16538_16540delATA, 16777T, 16823T, and 28736A) are novel mutations that did not appear in the Netherlands–Belgium or global datasets. The other 7 variants appeared but demonstrated low prevalence in both datasets (0.002%–0.130%). The most prominent of these was the 28139A mutation in a wastewater sample from March, which was detected in only 4 sequences worldwide and demonstrated both a strong temporal (all detected in March 2020) and regional signal (2 sequences from the Netherlands [EPI_ISL_422640 and EPI_ISL_422880], 1 from Denmark [EPI_ISL_444879], and 1 from Belgium [EPI_ISL_458209]).

18Finally, 4 variants (1440A, 11083T, 11109T, and 24862G) appeared at higher levels in both datasets (>0.5%); 11109T and 24862G were 28.5 and 14.3 times more prevalent in the Netherlands dataset than in the global dataset (Table 2). The other variants appeared at similar frequencies in all datasets.

19In addition to consensus sequences, LFV analysis is of value identifying potential local outbreaks. This identification could be achieved by detecting cluster-defining mutations that are associated with sequences from a particular geographic area. To associate the presence of a minor variant to sequences belonging to unique clusters, we mapped the 4 most prevalent LFVs onto the Netherlands–Belgium subsample and global subsample phylogenetic trees (Figure 3). For 3 variants (1440A, 11109T, and 24862G), the presence of the mutation and their clustering on the phylogenies were clearly associated. However, when 1 of these 3 variants was detected as an LFV in a sewage sample, the consensus sequence of this sample did not group with the cluster of clinical samples that contains the variant. For example, the 24862G variant in sample Tilburg-94339 was detected in 2 unique clusters within clade 20A, whereas its consensus sequence (hCoV-19/env/Netherlands/Tilburg-94339-I/2020) clustered within clade 20B, suggesting the presence of both clades in this sample. Although mutation 11083T was most prevalent in clade 19A, it was also scattered along the trees, indicating poor association with a particular clade.

20The use of wastewater sampling as a tool to learn more about the epidemiology and diversity of SARS-CoV-2 at a community level offers many advantages over human sampling. Sewage samples are relatively easy to collect, sampling bias toward severe cases does not occur, ethical issues are limited, and potentially fewer samples are required to determine temporal changes of viral infections in the community. Nevertheless, comprehensive comparisons with clinical surveillance are required to determine the extent and limits of using sewage as a surveillance or early-warning tool.

21We used nanopore and Illumina NGS analysis to study the diversity of SARS-CoV-2 in sewage and compared these results to the viral diversity found in clinical samples. To evaluate this diversity in a comprehensive fashion, we used the Nextstrain clade classification system because it is based on the use of signature mutations to assign sequences to a clade, enabling the association of SNPs or LFV to a particular clade, especially for genome sequences with <75% coverage.

22Our method enabled us to obtain complete or near-complete genomes from wastewater samples with Ct values of >5 Cts below the limit of detection and partial genomes for samples with higher Ct values. To increase the percentage of genome covered, a threshold of 10× coverage per position was used to generate consensus sequences from nanopore reads. The error rate with this threshold is <0.03%, and most of the mutations (132/145) listed have a coverage of >30×, which produces an error rate of 1/585,000 nt.

23Of note, we found sewage samples that clustered with sequences isolated from patients of the same region and LFV with a strong regional signal. In a recent study from the United States, wastewater contained SARS-CoV-2 genomes identical to those in clinical samples from the same region. Sewage samples can contain a mixture of SARS-CoV-2 viruses, which can be an indication of multiple viruses circulating within a community and perhaps in domestic and livestock animals. We applied a targeted amplification method and thus did not assess the presence of other viruses. Consensus sequence genomes from a wastewater sample can identify the predominant virus strain in a population, which is suitable for locations with few introductions of the virus. However, this approach is not appropriate for a population in which multiple virus strains are circulating in parallel. Moreover, it might lead to artificial consensus genomes that do not represent an existing virus.

24NGS analysis can unravel the diversity of viruses within a complex sample such as wastewater, particularly by using unbiased sequencing of the sewage virome. Nevertheless, the detection of variants of a virus in a single sample can be challenging because of the relatively low number of reads obtained for each virus. Targeted amplification and NGS of a small genome region of the virus of interest to determine the prevalence of virus variants within a single wastewater sample is more sensitive and less expensive; use of this approach has been reported for enteroviruses, human mastadenoviruses, and noroviruses. Because the diversity of SARS-CoV-2 is still limited, however, this approach would not be useful since no single small piece of the genome can reliably differentiate between clades or lineages. However, we demonstrated that some LFVs and SNPs can be linked to particular clusters or clades within trees without the need for a complete genome. To confidently determine the presence of a particular cluster within a sample, at least 2 LFVs associated with the cluster should be present at substantial levels. Furthermore, variant analysis can also be used to monitor the prevalence of biologically relevant mutations, such as D614G, which has been shown to increase infectivity in vitro and might be associated with higher transmission and death rates (; M. Cortey, unpub. data, https://www.biorxiv.org/content/10.1101/2020.05.16.099499v1). Within our dataset, clear temporal changes in the prevalence of LFVs or SNPs in sewage samples that correlated with changes in the clinical dataset were not detected during the first wave.

25The combination of whole-genome sequencing of clinical samples with epidemiologic data is vital for public health decision-making because it helps identify clusters of infection, new introductions of virus, and the expansion and decline of circulating strains. Cities with large numbers of visitors are expected to experience several introductions of the virus, whereas the opposite is expected for cities with low numbers of visitors. The use of NGS analysis of sewage samples to evaluate viral diversity within a geographic area and its changes over time can aid in decision-making. For example, in scenarios in which a large increase of viral diversity is detected in sewage, suggesting new introductions of virus, appropriate measures can be taken.

26Wastewater can also be used to monitor novel mutations. Our consensus and LFV analyses revealed 57 mutations that were not seen in the global database. These novel mutations might not have been detected for several reasons: they represent technical errors; the mutations did not stay within the population; or the mutations are associated with asymptomatic or mild disease, viruses from animal hosts, enteric shedding, or defective genomes. The presence of defective genomes has previously been suggested for the detection of LFVs that generate stop codons in clinical samples. Phenotypic studies could help determine the likelihood and biologic relevance of these novel mutations.

27In conclusion, this study illustrates the value of NGS analysis of wastewater to approximate the diversity of SARS-CoV-2 circulating in a community. Sequencing of wastewater samples could be a powerful tool to complement clinical surveillance or could be used independently in settings in which wide clinical sequencing is unfeasible. In addition, in-depth NGS analysis of wastewater samples can help in assessing changes in viral diversity, which can indicate the emergence of epidemiologically or clinically relevant mutations and thereby aid public health decision-making.