National Scale Real-Time Surveillance of SARS-CoV-2 Variants Dynamics by Wastewater Monitoring in Israel

In this report, we describe a national-scale monitoring of the SARS-CoV-2 (SC-2) variant dynamics in Israel, using multiple-time sampling of 13 wastewater treatment plants. We used a combination of inclusive and selective quantitative PCR assays that specifically identify variants A19/A20 or B.1.1.7 and tested each sample for the presence and relative viral RNA load of each variant. We show that between December 2020 and March 2021, a complete shift in the SC-2 variant circulation was observed, where the B.1.1.7 replaced the A19 in all examined test points. We further show that the normalized viral load (NVL) values and the average new cases per week reached a peak in January 2021 and then decreased gradually in almost all test points, in parallel with the progression of the national vaccination campaign, during February–March 2021. This study demonstrates the importance of monitoring SC-2 variant by using a combination of inclusive and selective PCR tests on a national scale through wastewater sampling, which is far more amendable for high-throughput monitoring compared with sequencing. This approach may be useful for real-time dynamics surveillance of current and future variants, such as the Omicron (BA.1, BA.2) and other variants.


Introduction
Since its emergence in December 2019, SARS -CoV-2 (SC-2) spread worldwide, causing the COVID-19 pandemic, with unprecedented impact on human lives, economy, social aspects, and public health throughout the world [1]. While SC-2 infection may cause a range of respiratory, gastrointestinal, and systemic symptoms, the majority of individuals show little to no symptoms [2][3][4][5][6][7] Clinical testing of individuals for SC-2 is the primary surveillance method for obtaining information for implementation of public health strategic interventions, such as quarantine, to mitigate the spread of the virus. The current gold standard for clinical testing is reverse transcriptase quantitative polymerase chain reaction (RT-qPCR) [3], which detects the viral RNA. Wastewater samples have been utilized to identify several pathogenic human viruses, and consequently, it has gained attention for assessing population-level trends of SC-2 infections.

Wastewater Sampling and Processing
The wastewater sampling sites (n = 13) selected for this study were collected from various regions of Israel covering more than 55% of Israel's population. These samples were collected from wastewater treatment plants (WWTP) between December 2020 and March 2021 (n = 117 samples). In each sampling cycle, 250 mL of untreated wastewater were collected every 30 min for 24 h by external automatic composite samplers. Samples were then transported to the laboratory under cooled conditions (4 • C). The sample material was concentrated as previously described [19], with minor changes. A duplicate homogeneous fraction of the wastewater sample (25 mL) was centrifuged (4696 g) for 5 min at 4 • C, and 20 mL of the supernatant was collected into 50 mL tube containing 0.26 M of MgCl 2 . The tubes were gently stirred for 5 min. Each sample duplicate was then passed through 0.45 µm pore-size, 47 mm diameter electronegative MCE membranes (Merck Millipore Ltd. Darmstadt, Germany) by peristaltic pumping. The membrane was transferred into a 50 mL tube containing 3 mL external lysis buffer (NucliSENS easyMAG, Marcy-l'Étoile, France). The tube with the membrane was gently stirred with the membrane inside for 30 min for virus inactivation and removal. The membrane was removed, and the lysis buffer containing the viral RNA was taken for RNA extraction.

Viral RNA Extraction
Total nucleic acid (NA) were extracted using the NucliSENS easyMAG system (Biomerieux, Marcy-l'Etoile, France) following the manufacturer's instructions. Briefly, 3 mL of lysis buffer containing the inactivated virus concentrated sewage were extracted using the easyMAG machine. Extracted NA were eluted using 55 µL elution buffer and stored at −70 • C.

Primers and Probes
Primers and probe for the OC43 spike control reaction were from Dare et al. [28]. Primers and probe for the SC-2 Envelope (E) gene reaction were based on the sequences published by Corman et al. [29], with the modifications detailed in Table 1. The primers for the SC-2 nucleocapsid (N) gene reaction were based on the CDC nCov-N1 reaction (https: //www.cdc.gov/coronavirus/2019-ncov/lab/rt-pcr-panel-primer-probes.html, (accessed 1 May 2022)), with the modification detailed in Table 1. The B.1.1.7 specific reaction was based on the CDC nCoV-N1, with modifications described previously [27]. The WT-specific N1 reaction was termed N WT , and the Alpha-specific reaction was termed N B.1.1.7 hereafter.

Design and Generation of In Vitro Transcribed Standard RNA Templates
In-vitro transcribed RNA molecules were used as control templates to maintain uniformity and to establish the reaction analytical limit of detection (LOD). Generation of RNA templates was performed using the Megascript kit (Thermo Fisher, https: //www.thermofisher.com/order/catalog/product/AMB13345#/AMB13345, accessed on 1 May 2022, Waltham, MA, USA), which allows in-vitro synthesis of RNA from DNA template using the T7 phage DNA-dependent RNA polymerase. Regions flanking the target sequences of each reaction were amplified using specific primers fused to the T7 phage minimal promoter sequence and were then used as template for the RNA synthesis. The resulting RNA templates were purified using the PSS magLEAD instrument (http://www.pss.co.jp/, accessed on 1 May 2022, Chiba, Japan). The RNA yield was determined using the NanoDrop spectrophotometer (https://www.thermofisher.com/i, accessed on 1 May 2022, Waltham, MA, USA) and used to generate standard curves for each reaction. The analytical limit of detection (LOD) was determined using serial dilutions of the RNA targets, with the multiplex mix. This calibration was then used to calculate the number of copies detected in wastewater samples.

Calculation of the Normalized Viral Load Value in Wastewater Samples
In order to compare the viral load as measured by qPCR, in different catchment points that have different flow rates and are draining regions with varying population size, we used a normalized viral load (NVL) conversion that takes into account these variables [30]. Using standard curves generated with the in vitro transcribed standard RNA templates, qPCR Cq values of the examined samples were converted into viral RNA copies. These values were then used for the NVL calculation according to the following equation: 1NVL (RNA copies population × day ) = RNA copies × Flow rate Population per 1000 where RNA copies per liter of wastewater is the converted from Cq value, flow rate is the measured flow rate (liter per day) as measured in the specific catchment point, and population per 1000 is the size of the population in thousands, which reside in the region that the specific catchment point is collecting. Calculation of the average NVL for each region was performed by determining the average NVL from each month (December, January, February, and March) for each WWTP. Then, the average NVL values from all the points of each region were summarized for each month. NVL values based the copy number derived from the E gene reaction were calculated to determine the total NVL of SC-2 in the sample pool. The percentage of the A19/B19 RNA and the B.1.1.7 RNA were determined by the NVL values of N WT reaction and N B.1.1.7 , respectively.

Whole Genome Sequencing and Bioinformatic Analysis
The procedure of library construction and WGS analysis was described in detail previously [27]. Briefly, the SC-2 genome was amplified using the Illumina COVID-seq kit, and libraries were prepared using the NexteraXT kit (https://emea.illumina.com/, accessed on 1 May 2022). Resulting Fastq sequence files that passed QC were compared to SC-2 reference genome (NC_045512.2) and analyzed for the presence and frequency of mutations that are associated with the B.1.1.7 variant. All genomic sequences used in this study were deposited in the GISAID database (https://www.gisaid.org/, accessed on 1 May 2022).

Ethical Statement
The sequencing of clinical samples was conducted according to the guidelines of the Declaration of Helsinki and was approved by the Institutional Review Board of the Sheba Medical Center (7045-20-SMC). Patient consent was waived as remains of clinical samples were used, and the data was anonymized during the analysis. The reaction that specifically detects the B.1.1.7 variant was described previously [27]. Briefly, a primer that binds specifically to the D3L substitution was designed and used, designated as 28257 N D3L VOC (highlighted in Table 1). This mutation is almost exclusive to variant B.1.1.7 and is not currently associated with any other major SC-2 variant, as evident from the NextStrain global mutation analysis ( Figure S1A). In order to perform a reciprocal reaction that detects only the WT lineage, a reciprocal primer that binds exclusively to the WT sequence was designed. The sequences and primer binding sites are illustrated in Figure S1B. In order to monitor the on-going spreading of variant B.1.1.7 in Israel, each sample was tested at least twice: once with the N B.1.1.7 assay and once with its reciprocal N WT reaction. For continuous monitoring of the variant circulation in wastewater, each sample was tested with a multiplex containing the E, N WT , and OC43 reactions and again with the E and N B.1.1.7 reactions. This enabled us to evaluate the overall SC-2 circulation in the examined wastewater treatment plants (WWTP) and, at the same time, determine the percentage of WT versus B.1.1.7 variants in that WWTP.

Development of RT-qPCR Assay for Detection of SC-2 in Environmental Samples
The designed primers were required to detect the presence of SC-2, to distinguish between the WT and variant B.1.1.7 strains, and to confirm the integrity of the sample preparation. In order to meet these requirements, the following triplex qPCR assay was developed. An inclusive SC-2 reaction targeting a conserved sequence within the SC-2 envelope (E) gene was used as a control for the presence or absence of SC-2 RNA in the sample. A second reaction targeting the N gene B.1.1.7 variant-associated mutation D3L was used as the differential reaction where the variant-specific primer was altered depending on the reaction specificity. A third reaction targeting the hCoV OC43 was used as an internal spiked control to evaluate the pre-extraction procedure integrity.
The specificity of the Alpha-specific reaction was already established previously [27]. In order to evaluate its selectivity with WW samples, 19A/19B and B.1.1.7 samples were examined each using either the N WT reaction or the N B.1.1.7 reaction. Figure S2 shows representative reactions in which the 19A/19B sample was positive for both the E and N targets when using the N WT reaction but only for the E target when using the N B.1.1.7 reaction. The B.1.1.7 sample gave the reciprocal result, i.e., positive for the E target only with the N WT reaction and positive for both targets with the N B.1.1.7 reaction ( Figure S2). The components of the three reactions were then combined together, and the sensitivity of each duplex was tested using in vitro synthesized RNA targets. Serial dilutions of the RNA targets were spiked into wastewater extraction to simulate an authentic sample background medium. The analytical limit of detection (LOD) derived from the dilutions testing was determined as follows: 35 copies per reaction for the E reaction, 12.9 copies for the N B.1.1.7 reaction, 13.5 copies for the N WT reaction, and 28 copies for the OC43 reaction ( Figure S3). For each test, the following positive controls were used: WT SC-2 RNA, B.1.1.7 RNA, and CoV OC43 RNA. SC-2-negative WW extraction was used as negative control. For quantitative analysis, serial dilutions of in vitro synthesized RNA targets were used.

Analysis of SC-2 Variants Dynamics Using the Differential Assay
Wastewater monitoring of SC-2 was performed from December 2020 to March 2021 in 13 WWTPs across Israel ( Figure 1). WWTPs population catchment area is estimated at about 5.2 million people, which is over 55% of Israel's population ( Table 2). The WWTPs were divided into four groups: north, central, Jerusalem, and south. In each regional group, 3-4 WWTPs were sampled periodically at least once a month. Each WWTP was sampled between 6 and 22 times depending on the sampling resources availability and the site accessibility. Each RNA extraction from each time point was tested at least twice. For each WWTP, the average normalized viral load (NVL) value was calculated from the repeats of the SC-2 E reaction. When applicable, the standard deviation was calculated and plotted on the representing graph. Analysis of the SC-2 dynamics in the northern region showed a mixed trend, where a slight decrease of the total normalized viral load (NVL) was observed in the El Hamra and Safed WWTPs, and a sharp decrease was evident in the Haifa WWTP. Turnover of the variant dominance from WT to B.1.1.7 was completed by the end of February 2021 ( Figure S4A). Monitoring of the four WWTPs of the central region showed a similar dynamic trend, with total NVL peaking in 21 January and slightly decreasing by March. As was in the northern region, by the beginning of March, nearly all samples contained only the variant B.1.1.7 ( Figure S4B). The total NVL values in the Jerusalem region WWTPs showed a similar course, in which the total NVL increased until January and then decreased by approximately 90% by March. The variant prevalence switched from 100% WT in December 2020 to 100% B.1.1.7 by the beginning of March 2021 ( Figure S5A). Interestingly, while the variant dynamics in the southern region were similar to those in the other three regions, the overall NVL showed a slight increase from December 2020 until the end of March 2021 ( Figure S5B). Comparison of the variant prevalence in February 2021 showed that the only WWTP in which the WT prevalence was above 46% was El Hamra (99%). In Safed, Ayalon, and Rahat, the B.1.1.7 variant consisted between 92% and 99% ( Figure S6).  In order to identify general trends in the variant dynamics, we clustered the results from the WWTPs of each region and examined the dynamics by region, as shown in Figure 2 and detailed in Table 2.  In order to identify general trends in the variant dynamics, we clustered the results from the WWTPs of each region and examined the dynamics by region, as shown in Figure 2 and detailed in Table 2. For each regional group, the NVL values from all WWTP sites within that group were used to calculate the regional average value. The average values obtained for each point were summarized for each region, and the percentage of the SC-2 variant was determined ( Figure 2). The variant turnover was similar in all four regions, switching from 100% WT in December 2020 to 100% variant B.1.1.7 by March 2021. The shift was somewhat delayed in the southern region, where in January, the WT prevalence was almost 100% compared to 15-45% in the other regions ( Figure 2). However, the shift in the southern region was faster compared to the Jerusalem region at 75% and 70% variant B.1.1.7 prevalence, respectively ( Figure 2). Although the switch from WT to B.1.1.7 was somewhat different in each region, all of them were 100% B.1.1.7 by mid-March (Figure 2, Figures S4 and S5). Comparison of The total NVL showed a moderate but clear decrease in the northern, central, and Jerusalem regions from January to March, while in the southern region, the NVL increased and remained constant during that period ( Figure 2). Detailed data for each WWTP site, including the calculated NVL values, are presented in Figures S4 and S5.

Variant Dynamics with Respect to Clinical Samples Sequencing and Vaccination Status
In order to assess the correlation between WW-based variant dynamics and se-  For each regional group, the NVL values from all WWTP sites within that group were used to calculate the regional average value. The average values obtained for each point were summarized for each region, and the percentage of the SC-2 variant was determined ( Figure 2). The variant turnover was similar in all four regions, switching from 100% WT in December 2020 to 100% variant B.1.1.7 by March 2021. The shift was somewhat delayed in the southern region, where in January, the WT prevalence was almost 100% compared to 15-45% in the other regions ( Figure 2). However, the shift in the southern region was faster compared to the Jerusalem region at 75% and 70% variant B.1.1.7 prevalence, respectively ( Figure 2). Although the switch from WT to B.1.1.7 was somewhat different in each region, all of them were 100% B.1.1.7 by mid-March (Figure 2, Figures S4 and S5). Comparison of The total NVL showed a moderate but clear decrease in the northern, central, and Jerusalem regions from January to March, while in the southern region, the NVL increased and remained constant during that period ( Figure 2). Detailed data for each WWTP site, including the calculated NVL values, are presented in Figures S4 and S5.

Variant Dynamics with Respect to Clinical Samples Sequencing and Vaccination Status
In order to assess the correlation between WW-based variant dynamics and sequenced clinical samples, the percentage of B.1.1.7-positive samples was calculated from the total number of samples that were sequenced and analyzed in the Israel Central Virology Laboratory (CVL) between December 2020 and March 2021 [27]. The results were then compared to the average percentages of B.1.1.7 and WT from all WWTP tested in the study. The average percent of either WT or B.1.1.7 was then plotted against the percentage of B.1.1.7-negative or B.1.1.7-positive clinical samples, respectively. As shown in Figure 3, on December 2020, 11% of the sequenced samples were identified as B.1.1.7, while only 1% of the total viral load calculated for WW samples was identified as such.  Figure 3). Notably, many of the clinical samples were not selected randomly for sequencing and may therefore somewhat skew the overall results. Nevertheless, these data indicated a clear correlation between the sequencing results and the dynamics inferred from WW examination.
During December 2020, when the Alpha variant was introduced into Israel, a national vaccination campaign was rolled out, using the two-dose BNT162b2 administration (Pfizer). By comparing the wastewater normalized viral load (WW NVL) of each region, we evaluated the possible effect of vaccination on the circulation of SC-2 in wastewater.
The total WW NVL were similar in the northern, central, and Jerusalem regions, with a NVL peak on January and a gradual decrease towards March, concurrent with the increase in second-dose vaccination rates. This trend was accompanied by a decrease in the new cases reported between January and March ( Figure 4 and Table 3). The average number of new cases in the northern, central, and southern regions was comparable, reaching approximately 1% in January 2021 and decreasing to 0.4-0.6% by March. The average number of new cases in the Jerusalem region was significantly higher, ranging from 3.4% in January to 0.8% in March (Table 3). Although the vaccination rates in the southern region were similar to the other regions, the number of new cases in that region increased, and the total WW NVL elevated accordingly (Figure 4 and Table 3). The decrease in total WW NVL and reported new cases in all other regions was parallel to the variant turnover, from WT to B.1.1.7 (Figures 2 and 4 Figure 3). Notably, many of the clinical samples were not selected randomly for sequencing and may therefore somewhat skew the overall results. Nevertheless, these data indicated a clear correlation between the sequencing results and the dynamics inferred from WW examination.
During December 2020, when the Alpha variant was introduced into Israel, a national vaccination campaign was rolled out, using the two-dose BNT162b2 administration (Pfizer). By comparing the wastewater normalized viral load (WW NVL) of each region, we evaluated the possible effect of vaccination on the circulation of SC-2 in wastewater.
The total WW NVL were similar in the northern, central, and Jerusalem regions, with a NVL peak on January and a gradual decrease towards March, concurrent with the increase in second-dose vaccination rates. This trend was accompanied by a decrease in the new cases reported between January and March ( Figure 4 and Table 3). The average number of new cases in the northern, central, and southern regions was comparable, reaching approximately 1% in January 2021 and decreasing to 0.4-0.6% by March. The average number of new cases in the Jerusalem region was significantly higher, ranging from 3.4% in January to 0.8% in March (Table 3). Although the vaccination rates in the southern region were similar to the other regions, the number of new cases in that region increased, and the total WW NVL elevated accordingly (Figure 4 and Table 3). The decrease in total WW NVL and reported new cases in all other regions was parallel to the variant turnover, from WT to B.1.1.7 (Figures 2 and 4).

Discussion
From the incursion of the SC-2 into Israel, until February 2021, the dominant SC-2 clade in Israel was the Wuhan 19B/20B. The Alpha variant was first identified in Israel at the end of December 2020 [27], and by March 2021, it consisted of more than 90% of randomly sequenced clinical samples (Israel Ministry of Health, data not published). Environmental surveillance proved as an important tool in monitoring the COVID-19 pandemic in populations and entire large geographic regions [21,31,32]. However, information on real-time dynamics of SC-2 variants on a national scale is currently insufficient. In this report, we describe the real-time prevalence of circulating SC-2 A19/B19 (Wuhan strain) and B1.1.7 variant, during a period of 4 months, across Israel. Some recent studies of SC-2 variants surveillance relied on the use of genomic sequencing or PCR detection of mutations that are not specific to a particular variant [33,34]. Other studies describe environmental surveillance using a combination of variant-specific and variant-nonspecific reactions [35,36]. Peterson et al. [36] describe the use of six reactions to indicate the presence of Alpha, Beta (B.1.351) or Gamma (P1) variant but not distinguish between them. The authors examined WW from five major cities as well as rural communities, determining the prevalence of the three mutations in each site. While the Alpha variant was verified by the specific N D3L-directed reaction, as we performed in our study, the other two reactions detected the spike 69-70 and N501Y mutations. These are characteristic but are not unique to the alpha variant [37]. Although the authors suggest that identifying these mutations may indicate the presence of Beta or Gamma variants, there was no definite evidence in the data presented that this was indeed the case. Furthermore, the need to run a separate assay for each mutation and the absence of an inclusive control (the "WT" reaction was run for each sample separately, in parallel with the "mutation" reaction) render this approach cumbersome. In our study, we combined an inclusive reaction (E-sarbeco) with the Alpha-specific reaction (D3L) and a control reaction (MS-2 phage spiking) to enable rapid and specific screening of a large number of samples.
The comparative analysis shows that the Alpha variant emerged in parallel on different sites in an apparently independent pattern. The first sites in which alpha variant was dominant by mid-February were remote from each other (Safed (Safed, Israel), Ayalon (Tel-Aviv, Israel, and Rahat (Rahat, Israel), Figure S4). This pattern further supports the assumption of independent Alpha variant spread. The different variant dynamics during January-February may reflect different infection events at the community level or simply result from sampling bias, as not all WWTPs were equally accessible for sampling. The complete turnover of all regions in March is consistent with the global variant dynamics, where the Alpha variant overtook the WT original SC-2 strain [37]. Comparison of the WW samples analysis to sequenced clinical samples showed an apparent preceding of the clinical samples turnover (Figure 3). This was likely due to the sequencing prioritization in December 2020 and January 2021, which included severe patients, incoming international travelers, samples with unusual PCR results, and other unique cases, all of which did not represent a random selection. Nevertheless, the overall variant prevalence, as reflected in clinical samples sequencing (WGS), was in agreement with the WW sampling results. In this respect, WW sampling may be more representative than individual sequencing since it is less prone to sampling bias, providing a pooled, heterogeneous sample material. These results indicated that WW monitoring could indeed be used as a good surveillance approach for identifying circulating viral variants.
Analysis of the variant dynamics with respect to population vaccination status showed a correlation between the increase in second-dose vaccination rates and the decline in total WW, NVL, and reported new clinical cases. This occurred regardless of the dominant variant. It is important to note that although the NVL dynamics in the southern region indicated an increase rather than a decrease, the calculated viral load was comparable to that of the northern and central regions and less than the Jerusalem region. The viral dynamics described in our study are consistent with the recent findings described by Yaniv et al. [35], which showed the turnover of the SC-2 variant during the beginning of 2021 in the city of Be'er Sheva. In a proceeding publication, the authors demonstrated that increased vaccination rates correlate with WW NVL decrease [37]. The authors conducted a thorough monitoring of the variant dynamics in Be'er Sheva and showed that the NVL increased from November 2020 to March 2021 and then started to decrease, following the increase in the second-dose vaccination rate [38]. We expanded this examination to 13 different sites across Israel and found a similar overall change with local differences between different sites, with evident correlation between the new cases and the NVL values.
In this study, which focused on the SC-2 WT-Alpha prevalence, we showed how a realtime and retrospective monitoring of WW could be used to establish the variant dynamics and reflect the population-level response to a national scale vaccination campaign over a period of 4 months.
Our data show a clear shift in the viral RNA circulation following the vaccination campaign, thereby demonstrating the importance of wastewater surveillance as an epidemiological tool. Previous reports were limited in the number of sites examined and mostly used sequencing to identify the circulation of variants of interest [39,40]. While sequencing provides complete information on the examined sample, it is expensive, not readily amendable for high throughput, and it is limited by the sample quality. We recently demonstrated that variant-specific PCR can determine SC-2 variant identity in cases of poor sample quality where sequencing fails [41]. Here, we show that variant-specific PCR can be a valuable tool to rapidly identify variant dynamics in WW in a cost-effective manner, which in turn enables a large-scale screening of WW samples to obtain important epidemiological information. Such an approach may be useful in current and future events of variant switch, as is happening now in many countries, with the BA.1 and BA.2 variants spreading (Bar-Or and Erster, manuscript in preparation). In addition to serve as an early warning system for the spreading of new variants, it allows a retrospective analysis of variant dynamics and vaccination response, thereby contributing to our understanding of the COVID19 progression and evolution.

Conflicts of Interest:
The authors declare no conflict of interests.

Ethical Statement:
The study was conducted according to the guidelines of the Declaration of Helsinki and was approved by the Institutional Review Board of the Sheba Medical Center (7045-20-SMC). Data regarding population vaccination and new positive SC-2 cases was obtained from the publicly available Israel Ministry of Health website at: https://datadashboard.health.gov.il/COVID-19/general (accessed on 1 May 2022).