Abstract
Wastewater-based epidemiology (WBE) is a recognised tool for tracking community transmission of COVID-19. From the second half of 2020, the emergence of new, highly infective, more pathogenic or vaccine-escape SARS-CoV-2 variants is the major public health concern. Variant analysis in sewage might assist the early detection of new mutations. Weekly raw sewage samples from 22 wastewater treatment plants (WWTPs) in Hungary (representing 40% of the population) were analysed between December 2020 and March 2021 for signature mutations N501Y and del H69/V70 of B.1.1.7 lineage by melting point genotyping and RT-digital droplet PCR (RT-ddPCR). The latter method proved to be more efficient in parallel detection of different variants and also provides quantitative information. Wastewater surveillance indicated that the B.1.1.7 variant first emerged in Budapest in early January 2021 and rapidly became dominant in the entire country. Results are in close agreement with the available clinical data (Pearson's correlation coefficient, R = 0.9153). RT-ddPCR was confirmed to be a reliable tool for tracking emerging variant ratios in wastewaters. It is a rapid and cost-effective method compared to whole-genome sequencing, but only applicable for the detection of known mutations. Efficient variant surveillance might require the combination of multiple methods.
HIGHLIGHTS
RT-digital droplet PCR was used to detect B.1.1.7 SARS-CoV-2 variant in sewage.
Signature mutations N501Y and del H69/V70 were detected in equal concentrations.
Signature mutations of B.1.1.7 became dominant within 6 weeks in Hungarian wastewater treatment plants.
Sewage data closely correlated with clinical emergence of B.1.1.7 cases.
Wastewater-based epidemiology is an efficient tool in tracing emerging SARS-CoV-2 variants of concern.
INTRODUCTION
The COVID-19 pandemic is a severe public health threat and its management is a serious challenge to governments. The need for rapid outbreak detection accelerated the introduction of new methods. Wastewater-based epidemiology (WBE) has been used for decades to monitor the health of the population (Sims & Kasprzyk-Hordern 2020), but it gained global attention in the current COVID-19 pandemic as an important complementary method to clinical surveillance.
Although it is a respiratory virus, SARS-CoV-2 is shed in faeces and (rarely) in urine by both symptomatic and asymptomatic individuals. Shedding may start before clinical symptoms and can extend up to several weeks, especially in faeces (Walsh et al. 2020). The presence of SARS-CoV-2 in the sewage and the correlation of its concentration to COVID-19 incidence rates were investigated in many countries by several research groups (e.g. Medema et al. 2020; Weidhaas et al. 2021; Wu et al. 2020). Most studies find a correlation between the SARS-CoV-2 concentration in wastewater and the number of clinically confirmed new COVID-19 cases with an offset of 4–10 days after wastewater sampling. The role of WBE in the COVID-19 pandemic and in possible subsequent epidemics is acknowledged by the WHO (2020). The European Commission issued a recommendation to its Member States to introduce wastewater surveillance for SARS-CoV-2, covering at least the agglomerations over 150,000 inhabitants (European Commission 2021). The recommendation also calls for the determination of virus variants from sewage.
The emergence of new SARS-CoV-2 variants is the main driver of consecutive new waves of the pandemic. The variants that pose an increased risk to global public health are termed variants of interest (VOIs) and variants of concern (VOCs) by the WHO (2021). The WHO in February 2021 defined VOIs as genetic changes that are likely to affect virus characteristics, such as transmissibility, disease severity, immune or treatment escape, and its epidemiological impacts suggest an emerging risk to global public health. VOCs are VOIs that have been demonstrated to be associated with significantly increased transmissibility, virulence or decreased effectiveness of diagnostics, vaccines or therapeutics. According to the WHO, the threshold for the determination of VOCs is high in order to focus attention and resources on the variants with the highest public health implications.
B.1.1.7 lineage or Alpha variant was the first VOC emerging in Europe from December 2020 (Davies et al. 2021a). The first Hungarian case was confirmed on 13 January. The lineage was defined with 17 non-synonymous mutations. The most important ones include N501Y (related to increased ACE2 receptor binding affinity) and deletion 69H/70 V in the S gene (Rambaut et al. 2020) encoding the surface spike protein of the virus. Increased transmissibility and increased disease severity were both demonstrated for this lineage (Davies et al. 2021b).
While lineage B.1.1.7 became dominant by the beginning of February 2021 in most European countries, B.1.351 and P.1 variants also appeared in the continent, but did not reach high prevalence in the first half of 2021 (ECDC 2021). These lineages also contain the N501Y amino acid change, but no 69H/70 V deletion. Both were declared VOC as they were considered vaccine-escape variants since they harbour the E484 K mutation (Wang et al. 2021a, 2021b).
Estimating the proportion of variants among the COVID-19 cases requires continuous and representative clinical investigation, which can be very resource-intensive. WBE was suggested as a cost-effective complementary method for the surveillance of virus variants. Whole-genome sequencing (WGS) is the generally accepted method of variant identification. Jahn et al. (2021) could detect 10 out of the 17 B.1.1.7 mutations with this method in a ski resort in Switzerland, 1 day before the first reported human case. Other methods are also proposed for variant detection from wastewaters. Allele-specific PCR methods screen for one or few signature mutations. Heijnen et al. (2021) used digital droplet PCR for selective quantification for N501Y, a key mutation in B.1.1.7, B.1.351 and P.1. Allele-specific RT-qPCR (Lee et al. 2021) and nested RT-PCR (La Rosa et al. 2021) were also suggested to be applicable for this purpose.
The aim of the present study is to identify a suitable method to assess the prevalence and proportion of SARS-CoV-2 VOC lineages in Hungary. Hungarian wastewater samples from December 2020 to mid-March 2021 – the period when the lineage B.1.1.7 variant became dominant in Hungary – were screened by multiple methods. A method applicable for the rapid detection of variants from sewage will be an important asset in the preparation for further pandemic waves.
MATERIALS AND METHODS
Sample collection
Raw sewage samples were collected in the wastewater treatment plants (WWTPs) of Budapest (n = 3), of all Hungarian county seats (n = 18), and a composite sample from five cities in agglomeration of Budapest (Tököl, Biatorbágy, Szigetszentmiklós, Budakeszi and Százhalombatta). Six sampled WWTPs (Budapest Central, Kaposvár, Miskolc, Nyíregyháza, Székesfehérvár and Szombathely) were equipped with built-in composite auto-sampler; from the other sampling points, grab samples were taken at the peak load of the WWTP in the morning hours. Samples were collected by accredited personnel of the WWTPs and transported within 24 h at 4 °C to the laboratory where it was processed within 6 h. Weekly samples from 1 December 2020 (Budapest) or from 15 January 2021 (other sampling points) to 15 March 2021 were investigated in this study (51 composite and 158 grab samples in total). The population served by each WWTP is included in Supplementary Material, Table S1.
Sample preparation and concentration
Sewage samples were concentrated by an in-house ultrafiltration method as described previously (Róka et al. 2021). Briefly, samples were sedimented by centrifugation (4,500 g, 30 min, 4 °C), followed by membrane ultrafiltration (Suez Water Technologies & Solutions, Membrane Research Center, Tatabánya, Hungary). In this process, 50 ml sewage is concentrated to 1 ml. External process control (wastewater sample spiked with inactivated SARS-CoV-2) was used to ensure consistency (Róka et al. 2021).
RNA extraction
For melting curve genotyping and RT-digital droplet PCR (RT-ddPCR), RNA was isolated using a commercial kit (QIAamp Viral RNA Mini Kit, Qiagen, Germany) in accordance with the manufacturer's instructions. The method yields 30 μl RNA from 140 μl concentrate.
Detection of SARS-CoV-2
For the detection of SARS-CoV-2 RNA, quantitative RT-PCR was conducted on the LightCycler 480 Instrument II (Roche, Switzerland) platform using the LightCycler Multiplex RNA Virus Master kit (Roche, Switzerland). Primers and probes used were specific for the nucleocapsid 1 (N1) gene of the virus (2019-nCoV-N1-F: GAC CCC AAA ATC AGC GAA AT; 2019-nCoV-N1-R: TCT GGT TAC TGC CAG TTG AAT CTG; 2019-nCoV-N1-P: FAM-ACC CCG CAT TAC GTT TGG TGG ACC-BHQ1) (CDC, Atlanta, https://www.cdc.gov/coronavirus/2019-ncov/lab/rt-pcr-panel-primer-probes.html). Virus copy number was calculated based on the standard curve of a commercially available copy control (EURM-019) obtained from the European Commission Joint Research Center, Geel, measured with an automatic second derivative maximum method threshold. The limit of detection (LOD) of the RT-PCR was measured using log10 serial dilution of the EURM-19 copy control.
Melting curve genotyping
For qualitative detection of spike mutations N501Y and del H69/V70, the VirSNiP SARS-CoV-2 Spike N501Y and Spike del H69/V70 assays were used (TIB MOLBIOL, Berlin, Germany). The genotyping RT-PCR and melting analysis were conducted on the LightCycler 480 Instrument II platform using LightCycler Multiplex RNA Virus Master kit in 20 μL final reaction volume and using simple probe detection format. For the melting curve analysis, Tm calling was performed using the Standard ModularDX Program for Roche LC480 Instruments. The Tm temperature of the samples was compared to the N501/501Y and del H69/V70 positive and negative controls and to the RT-PCR negative controls.
RT-ddPCR for mutation detection
For simultaneous quantitative detection of the signature mutations of the B.1.1.7 variant and wild-type of the same positions, ddRT-PCR mutation detection assays were used (Bio-Rad, CA, USA). For 501Y and the wild-type variant (N501) (assay ID: dMDS731762551), for the del69H/70 V and the wild-type variant (no deletion) (assay ID: dMDS284738817), primers-probes were used. The 22 μl reaction volume from One-Step RT-ddPCR Advanced Kit for Probes contained 5.5 μl one-step RT-ddPCR Supermix, 2.2 μl reverse transcriptase, 1.1 μl 300 mM DTT, 1 μl mutation detection assay, 7.2 μl water and 5.0 μl purified RNA. Droplets were generated using the Automated Droplet Generator (Bio-Rad, CA, USA), and RT-ddPCR was carried out on the Bio-Rad T100 Thermal Cycler platform (25 °C 30 min pre-incubation, 50 °C 60 min reverse transcription, 95 °C 10 min enzyme activation, 40 cycles of 95 °C 30 s denaturation and 55 °C 60 s annealing/extension; 98 °C 10 min enzyme deactivation and overnight droplet stabilisation at 4 °C; ramp rate was decreased to 2.0 °C/s in all steps). The RT-ddPCR product read-out was performed using the QX200 system (Bio-Rad) and analysed using the QuantaSoft™ Analysis Pro 1.0.596 software (Bio-Rad). Wild-type, B.1.1.7 lineage and non-template control were used in every reaction.
The LOD and the limit of quantification (LOQ) were determined using replicates of 10-fold dilutions of purified SARS-CoV-2 RNA (same applied as positive amplification control in RT-qPCR). The LOD was defined as the last serial dilution detected in 95% of replicates and LOQ as the lowest dilution showing a coefficient of variation percentage below the threshold (CV% = 25) (Rocchigiani et al. 2020). The calculated value both for the LOD and the LOQ was 0.703 copies/μl (Supplementary Material, Table S3).
Sewage data
The Enterococcus count of the sewage samples was determined according to ISO 7899-1 standard methods. MUD/SF plates were obtained from Bio-Rad (CA, USA). Plates were incubated at 44 °C for 48 h and read under UV light.
Sewage flow data (daily volume) were provided by the WWTP operators (data shown in Supplementary Material, Table S2).
Epidemiological data
The total number of new COVID-19 cases was obtained from the National Surveillance Database of National Public Health Centre, and the number of VOC-infected cases (B.1.1.7, B.1.351 and P.1) was obtained from the database of the COVID Laboratory Working Group of the National Public Health Centre.
Statistical analysis
Microsoft Excel 2010 (version: 14.0.7232.5000) was used for the statistical analysis of data. Pearson's correlation coefficient was calculated for comparing the concentration of the two investigated mutations and investigating the association of clinical B.1.1.7 cases and the variant's concentration in sewage. Linear regression was used to confirm the association between the SARS-CoV-2 gene copy numbers in sewage and the weekly number of new cases.
RESULTS
SARS-CoV-2 concentration
During the monitored period, the SARS-CoV-2 concentration in the wastewater samples decreased slowly until the third week of 2021, parallel to the decline of the second wave. From week 4 in 2021, a rapid increase was observed (Figure 1; Supplementary Material, Table S2). The comparison of the SARS-CoV-2 genome copy numbers and the weekly average of new cases by linear regression confirmed that trends sewage data precede the changes in the case number by approximately 1 week (R = 0.843, p = 0.0002).
Population weighted average of SARS-CoV-2 concentrations and 7-day rolling average (from Thursday to Wednesday) of new COVID-19 cases.
Population weighted average of SARS-CoV-2 concentrations and 7-day rolling average (from Thursday to Wednesday) of new COVID-19 cases.
Melting curve genotyping
Melting curve genotyping is only suitable for the qualitative detection of mutations. However, if different variants are present, multiple melting temperatures are also detectable. Wastewater samples spiked with heat-inactivated SARS-CoV-2 (using confirmed B.1.1.7 lineage and wild-type viruses in a final concentration of approximately 106 GC/l) were used to confirm the suitability of the method for wastewaters. During melting curve analysis, multiple peaks of melting temperatures were observed in sewage samples spiked with both isolates (B.1.1.7 and wild-type in 1:1 proportion).
Twenty-one raw sewage samples (6 composite and 15 grab samples) from the 6–7th week of 2021 were analysed with this method. Wild-type N501 was detected in two samples, N501Y genotype in nine samples and both variants in two samples. In eight cases, N501Y detection PCR was unsuccessful. Screening for del H69/V70 was only successful in five cases, of which two samples showed wild-type genotype, two samples contained the deletion and one sample contained both variants (Supplementary Material, Table S2).
RT-digital droplet PCR
The method allows parallel quantification of the mutant and wild-type of a genomic locus in the S gene (N501Y or delH69/V70). 212 samples (47 composite and 165 grab samples) obtained between weeks 49/2020 and 10/2021 (in the case of Budapest) and between weeks 3/2021 and 10/2021 (at other sampling sites) were analysed by RT-ddPCR.
The LOD for different variants – as a percentage of the total concentration – was estimated based on the process control samples. To determine the virus concentration of the process control in the wastewater, a natural wastewater sample spiked with heat-inactivated, wild-type SARS-CoV-2 was used. By using this method, we were able to efficiently recover other lineages present naturally in low concentrations (0.5–5%) in the spiked sample.
Wastewater samples from Budapest were investigated for a longer period, as earlier emergence of new types was expected in the capital based on clinical data. The results show high levels of the wild-type virus at the beginning of December. The concentration of the wild-type strain started to decrease from mid-December and continued to decrease over the study period. The first part (until week 4 of 2021) of the observed trend is associated with the decline of the second wave of the outbreak but after that, the cumulated SARS-CoV-2 concentration started to increase (see Figure 1). Hence, the further decrease of the wild-type virus is due to the emergence of other variants. Both of the monitored spike mutations occurred at low (but not zero) concentrations at the beginning of December and increased rapidly from the beginning of February (Figure 2(a)). Comparing the proportion of the mutant and wild-type variants, the first emergence of N501Y and del H69/V70 can be observed from the third week of 2021 (Figure 2(b)), and by the end of February, it became the variant almost exclusively present in wastewaters (Figure 2(b)).
(a) Weekly changes of the emergence of signature mutations of the S gene specific for the B.1.1.7 variant in Budapest (average if three WWTPs). (b) Changes in the ratio of the wild-type and presumptive B.1.1.7 SARS-CoV-2 variants in Budapest wastewater (average of three WWTPs) between December 2020 and March 2021. Presumptive B.1.1.7 genome copy numbers were calculated as the average of N501Y and del H69/V70 copy numbers, while wild-type was expressed as the average of the non-mutated copy numbers of the two sites.
(a) Weekly changes of the emergence of signature mutations of the S gene specific for the B.1.1.7 variant in Budapest (average if three WWTPs). (b) Changes in the ratio of the wild-type and presumptive B.1.1.7 SARS-CoV-2 variants in Budapest wastewater (average of three WWTPs) between December 2020 and March 2021. Presumptive B.1.1.7 genome copy numbers were calculated as the average of N501Y and del H69/V70 copy numbers, while wild-type was expressed as the average of the non-mutated copy numbers of the two sites.
Similar trends were observed in every WWTP sampled in the national surveillance (Table 1). The average proportion of variants N501Y and del H69/V70 was nearly zero in mid-January and become almost exclusively present by mid-March. However, several weeks of discrepancy are observed in the spread of the new variant across the country.
Average of N501Y and del H69/V70 proportion in the investigated WWTPs
![]() |
![]() |
□ no data. The proportion of presumptive B.1.1.7 lineage in sewage sample (average ratio of N501 and delH69/V70).
*The mixed sample from five cities (Tököl, Biatorbágy, Szigetszentmiklós, Budakeszi and Százhalombatta).
The cells are coloured on a colour scale indicating the presence of different variants (green, only wild type is present; from light green to orange, both types are present, with increasing ratio of mutants; red, only mutants are present). Please refer to the online version of this paper to see this table in colour: http://dx.doi.10.2166/wh.2022.179.
B.1.1.7 cases and signature mutations in the wastewater
The ratio of N501Y+del69H/V70 mutations in the clinical samples analysed by the National Public Health Center COVID Laboratory Working Group was similar to that observed in wastewater in the investigated period (Figure 3). Clinical samples positive for N501Y+del69H/V70 mutations by the VirSNiP Mutation Assays were also confirmed partially by WGS, indicating the B.1.1.7 variant present in 100% of the tested samples. Presumptive lineage B.1.1.7 was accountable for approximately 10% of all cases in mid-January 2021 and reached 90% by mid-March.
Percentage of presumptive lineage B.1.1.7 (based on N501Y and del H69/V70 mutations screening) in clinical samples and the same mutations measured in wastewater (national average).
Percentage of presumptive lineage B.1.1.7 (based on N501Y and del H69/V70 mutations screening) in clinical samples and the same mutations measured in wastewater (national average).
DISCUSSION
In the present study, two methods were tested for the detection of SARS-CoV-2 VOCs in wastewater, melting curve genotyping and RT-ddPCR. The detection of N501Y or wild-type N501 was successful with both methods from 13 samples. The two samples showing wild-type N501 according to melting curve genotyping contained 6–13% N501Y according to RT-ddPCR, while N501Y positive samples contained 56–93% N501Y. The two samples with multiple maximums of melting temperature showed 64 and 94% N501Y according to RT-ddPCR. Similar results can be observed by the H69/V70 deletion screening methods (Supplementary Material, Table S2).
The advantages of the melting curve genotyping method are the low cost, rapid identification and easy access to the materials. For these reasons, investigations were carried out nearly in real-time in February, from fresh samples. However, low sensitivity and lack of quantitative data are strong limitations of melting curve genotyping.
The results of RT-ddPCR indicate a rapid dispersion of the two signature mutations of the B.1.1.7 variant in the country. The variants first emerged in Budapest and its agglomeration but appeared in all WWTPs within 4 weeks. The onset of the emergence of the variants did not show a clear correlation either with the size of the population serviced by the WWTP or the geographical location. The increasing trend was consistent in most WWTPs, though there were occasional outliers. The reason behind unexpectedly high (e.g. Szeged, 3rd week or Győr, 4th week) or unexpectedly low (e.g. Debrecen, 9th week, Veszprém, 8th week or Zalaegerszeg, 10th week) VOC ratios is not clear. One potential explanation is the inherent heterogeneity of the wastewater which may lead to a bias, especially in the case of the grab sampling. Outlying results were more frequent in those WWTPs where an automated sampler was not available.
The frequency of B.1.1.7 variant cases and the ratio of its signature mutations in sewage increased rapidly between mid-January and mid-March 2021 (see Figure 3). Though clinical results are not representative for the entire country (figures are based on clinical samples tested in NPHC, which represents approximately 3–15% of the new daily cases in the country), the correlation is statistically significant (Pearson's correlation coefficient, R = 0.921, p<0.05). Other studies using RT-ddPCR (Heijnen et al. 2021) or sequencing (Jahn et al. 2021) for variant analysis in sewage have also indicated that the ratio of different SARS-CoV-2 variants represents well the distribution of different lineages in clinical cases.
The two investigated mutations are key changes in the B.1.1.7 lineage but it is not defined exclusively by them: N501Y can also be found in B.1.351 and P1 lineages, while B.1.1.7 also harbours several other non-synonymous mutations (Rambaut et al. 2020). The RT-ddPCR method detects the two mutation sites in two separate reactions; therefore, it is not automatically evident that they are in present the same viral genomes. However, the prevalence of N501Y and del 69H/70 V was closely linked in the investigated samples (Pearson's correlation, R = 0.9153, p<0.05; Figure 4).
Correlation between N501Y and del H69/V70 concentrations (GC/μl, in RT-ddPCR reaction).
Correlation between N501Y and del H69/V70 concentrations (GC/μl, in RT-ddPCR reaction).
A less clear correlation was observed between the quantitative RT-qPCR (CDC N1 primer-probe set) and RT-ddPCR (S gene) results (Pearson's correlation, R = 0.7173, p<0.05; Figure 5). Heijnen et al. (2021) found a very similar correlation between the two datasets (RT-qPCR with N2 primer set compared to RT-ddPCR with N501Y/N501 WT primer set). The total concentration of different variants measured by RT-ddPCR was lower than the qPCR results by an average of 62% (average of N501Y+N501 WT and del H69/V70+wild-type without deletion). The difference is probably due to the different quantification methods and the different target genes of the two types of PCR. A similar difference was seen in RT-ddPCR and RT-qPCR comparison with the same primer set (N1, data not shown).
Correlation between RT-qPCR results (N1 primer set) and RT-ddPCR (average of N501Y+N501WT and del H69/V70+WT without deletion).
Correlation between RT-qPCR results (N1 primer set) and RT-ddPCR (average of N501Y+N501WT and del H69/V70+WT without deletion).
CONCLUSION
Results of the present study confirmed that RT-ddPCR-based methods are rapid and reliable tools for tracking VOCs in sewage, and the results are also indicative for community circulation of different variants. It is more convenient for surveillance than sequencing due to its lower time demand (6 h vs. several days) and cost-efficiency, especially in resource-limited settings.
Though RT-ddPCR uses signature mutations for variant detection, and not the full mutation profile, the close correlation of the two mutations' concentration and the association with the clinical data strongly suggest that their presence was due to the emergence of the B.1.1.7 variant. WGS and quasi-species analysis could provide further confirmation.
The tested melting point genotyping method has limited potential for the analysis of variants in sewage and should only be used in the absence of the other available options.
Recent developments in PCR-based methods for VOC screening in wastewaters offer solutions for the limitation of the currently available techniques. Lee et al. (2021) recently published a method for allele-specific RT-qPCR for three genomic loci predictive for B.1.1.7 (HV69/70del, Y144del and A570D), while La Rosa et al. (2021) used nested RT-PCR assays for rapid variant screening. These recently published methods seem to be more sensitive for mutations, and the RT-qPCR method is also suitable for quantitative analysis. However, rapid, reliable and quantitative variant analysis in sewage might still require a combination of different techniques (e.g. allele-specific RT-qPCR or RT-ddPCR and WGS).
DATA AVAILABILITY STATEMENT
All relevant data are included in the paper or its Supplementary Information.