Introduction

SARS-CoV-2, the etiological agent responsible for the deadly respiratory disease known as COVID-19, emerged in Wuhan (China) in early December 2019. It was on March 11, 2020, when the World Health Organization (WHO) declared the novel coronavirus outbreak a global pandemic (WHO 2020). The transmission dynamics of COVID-19 is due to two mechanisms: human-to-human and air pollution-to-person transmission. The first one is through respiratory droplets and direct contact with infected people, and it is influenced by the population density. However, Coccia (2020) highlighted the importance of the air pollution in the rapid spread of COVID-19. They suggested that climate and meteorology factors such as pollution or wind which affect the content of microbes in the atmosphere, accelerating the spread of SARS-CoV-2 in the environment, including wastewater. Besides, although the fecal–oral transmission of SARS-CoV-2 has not been proved yet (Albert et al. 2021; Guo et al. 2021), COVID-19 dissemination could increase due to inadequate sanitation facilities or poor fecal sludge management in hospitals environments (Amin et al. 2023).

The viral surface spike glycoprotein S plays an important role in the infection process, mediating viral entry into host cells through the union of the spike receptor-binding domain (RBD) and the angiotensin-converting enzyme2 (ACE2) receptor (Lan et al. 2020). ACE2 is abundant in glandular cells of the gastric, duodenal, and rectal epithelium (Xiao et al. 2020), where the virus persists longer than in the respiratory tract (Hu et al. 2020; Zheng et al. 2020). Therefore, people infected with SARS-CoV-2 shed the virus in their feces despite being asymptomatic (Lai et al. 2020; Rothe et al. 2020) or having tested negative in nasopharyngeal samples (Jiang et al. 2020; Wu et al. 2020). It has been shown that the virus can remain in stool samples for up to 5 weeks after the onset of symptoms (Wu et al. 2020). As a result, SARS-CoV-2 RNA from most infected patients ends up in wastewater treatment plants (WWTPs). Therefore, the analysis of wastewater samples can provide effective epidemiological surveillance.

Wastewater-based epidemiology (WBE) has been used as a method of early detection and direct mitigation of poliovirus outbreaks in Israel and Egypt (Blomqvist et al. 2012; Brouwer et al. 2018; Kopel et al. 2014) or Norovirus and Hepatitis A in Sweden (Berchenko et al. 2017; Duintjer Tebbens et al. 2017; Hellmér et al. 2014). Currently, RT-qPCR assays are used for clinical trials and the detection of viral RNA in WWT (Acosta et al. 2022; Kolarević et al. 2022; Tanimoto et al. 2022; Yanaç et al. 2022).

In most WBE SARS-CoV-2 studies, a significant correlation has been found between the viral load measured in wastewater and clinical cases of COVID-19 (Maida et al. 2022; Pillay et al. 2021; Vallejo et al. 2022), which demonstrates that an increase on viral load in WWT can warn about the emergence of a new variant or an outbreak (Barua et al. 2022; Daleiden et al. 2022; Kuhn et al. 2022; Li et al. 2022; Monteiro et al. 2022; Padilla-Reyes et al. 2022; Robotto et al. 2022; Sangsanont et al. 2022; Wu et al. 2022). In addition, WBE has been used to estimate the number of people infected with SARS-CoV-2 in the population (Chavarria-Miro et al. 2021; McMahan et al. 2021; Saththasivam et al. 2021; Tharak et al. 2022; Vallejo et al. 2022). However, there are few studies focused on monitoring the virus and predicting COVID-19 clinical cases in specifically rural areas or at the building level (Jarvie et al. 2023). It should be noted that the effectiveness of the post-factum methods in WBE studies, that is, an early warning analysis based on previously reported both wastewater and clinical data, has been proved (Zhao et al. 2023). The implementation of both post-factum and real-time early warning methods in a parallel and complementary way could increase the precision of the WBE analysis by almost 100%.

The emergence of new SARS-CoV-2 variants since late 2020 prompted the classification by the European Centre for Disease Prevention and Control (ECDC) of Variants of Interest (VOI) Variants of Concern (VOC), Variants Under Monitoring (VUM), and de-escalated variants (ECDC 2020). VOCs have a significant epidemiological impact due to their increase in transmissibility or in the severity of the disease they cause (Gobeil et al. 2021; Harvey et al. 2021; Riou et al. 2022). From the beginning of the pandemic up to now, Alpha (B.1.1.7), Beta (B.1.351), Gamma (P.1), Delta (B.1.617.2), and Omicron (B.1.1.529) VOCs have been most relevant causing different epidemic waves. The emergence of the variants has made the need for sequencing in wastewater more evident. On March 17, 2021, the European Commission (EC) published several recommendations to foster the surveillance of SARS-CoV-2 and its variants in wastewater (EC 2021a). Not surprisingly, genomic surveillance and monitoring are highly effective in detecting the emergence of SARS-CoV-2 variants in imported cases and in small outbreaks (Fontenele et al. 2021; Pechlivanis et al. 2022; Pérez-Cataluña et al. 2022; Peterson et al. 2022; Sapoval et al. 2021). Moreover, our team has been previously focused on estimating the proportions of variants in the population based on mutations data found in wastewater samples using statistical models (López de Ullibarri et al. 2023). Similar models have also been described by other authors (Gafurov et al. 2022; Pipes et al. 2022; Valieris et al. 2022).

The COVIDBENS project was one of the earliest WBE on COVID-19 projects in the world, starting on April 14, 2020 (Vallejo et al. 2022). In the present work, we analysed 863 wastewater samples from June 2020 to March 2022 collected from the WWTP Bens (A Coruña, Spain), a public company that serves a population of ca. 400.000 inhabitants. The main objectives of this multidisciplinary study were the following: (1) to monitor at real time the COVID-19 pandemic by tracking SARS-CoV-2 viral load in wastewater in the metropolitan area of A Coruña; (2) to estimate the number of people infected with SARS-CoV-2 in the local population, including symptomatic and asymptomatic people, using statistical models described before by our team (Vallejo et al. 2022); (3) to act as an early warning system for predicting new outbreaks before the health system; (4) to detect SARS-CoV-2 mutations in wastewater by next-generation sequencing (NGS); (5) to estimate the percentage of SARS-CoV-2 variants based on the mutations found in the wastewater using previous statistical models described by our team (López de Ullibarri et al. 2023); and (6) to inform the authorities and all citizens about the real evolution of the COVID-19 pandemic providing a public service. In this manuscript, we describe how COVIDBENS alerted about the increase of SARS-CoV-2 in the population on multiple occasions, predicting the most significant outbreaks during the COVID-19 pandemic and the emergence of the different variants that affected this area. COVIDBENS has been extremely useful as an early warning system for decision-making and for evaluating the impact of SARS-CoV-2 control measures or vaccination campaigns in our local area, A Coruña (Spain).

Materials and methods

Sample and data

Wastewater sampling

The WWTP of Bens serves a population of ca. 400.000 inhabitants from the metropolitan area of A Coruña (NW, Spain) that includes the municipalities of A Coruña, Arteixo, Cambre, Culleredo, and Oleiros (Fig. 1). Composite sewage samples were collected twice or thrice a week using automatic samplers installed by operators both at the entrance (influent) of the WWTP and at the sewers of each municipality. These automatic samplers were programmed to collect 150 mL of wastewater every 15 min over a 24-h period, resulting in a 600 mL bottle per hour. Then, the resulting 24 bottles were merged into a larger one, and a sample of 100 mL was kept on ice and processed immediately after reception at the lab. For the present study, samples from the WWTP were collected from June 4, 2020, to March 17, 2022, and, in the case of A Coruña, Arteixo, Cambre, Culleredo, and Oleiros, were collected separately from December 24, 2020, to March 17, 2022.

Fig. 1
figure 1

Map showing the sampling area served by the WWTP of Bens (▲), including the areas of the municipalities of A Coruña, Arteixo, Cambre, Oleiros, and Culleredo which correspond to individual sewers

Full size image

Sample processing

The 100 mL samples were centrifuged for 30 min at 4000 × g at 4 °C and filtered through 0.22-μm pore membranes (Merck Millipore, USA). Pellets were discarded, and the supernatants were concentrated and dialyzed with SM buffer (50 mM Tris–HCL, 100 mM NaCl, and 8 mM MgSO4) using the Vivaspin Turbo 15 with a polyethersulfone membrane of 30-kDa (Sartorius, Germany). Finally, samples were preserved in 500 μL of RNAlater™ solution (Thermo Fisher Scientific, USA) at − 80 °C for further analyses.

For viral load determination, SARS-CoV-2 RNA was isolated from 100 μL of the thawed samples using the QIAamp Viral RNA Mini Kit (Qiagen, Germany), according to the manufacturer’s protocol. Samples were first mixed with AVL buffer supplemented with carrier RNA to ensure the binding of viral RNA to the QIAamp silica-based membrane. Then, samples were washed twice to remove PCR inhibitors and eluted in 70 μL of RNase-free water (Thermo Fisher Scientific, USA). For sequencing, RNA was extracted again from the remaining amount of sample (400 μL) and eluted in 80 μL of RNase-free water (Thermo Fisher Scientific, USA). Extracted RNA was stored at − 80 °C until use.

SARS-CoV-2 viral RNA was detected by one-step real-time reverse transcription quantitative PCR (RT-qPCR) using the TaqPath COVID-19 RT-PCR Kit (Thermo Fisher Scientific, USA) and a CFX96 Thermal cycler (Bio-Rad, USA). This kit includes primer pairs targeting SARS-CoV-2 specific genome regions such as ORF1ab, the N gene, and the S gene. Recommendations for minimising errors in RT–PCR detection were followed (Ahmed et al. 2022). Reactions were initiated with an uracil N-glycosylase (UNG) incubation at 25 °C for 2 min and a reverse transcription at 53 °C for 10 min, followed by a polymerase activation at 95 °C for 2 min; 40 cycles of 95 °C for 3 s for DNA denaturation; and 60 °C for 30 s for annealing and extension. A 25 μL reaction contained 15 μL of reaction mix (6.25 μL of TaqPath 1-Step Multiplex Master Mix, 1.25 μL of COVID-19 Real-Time PCR Assay Multiplex, and 7.5 μL of nuclease-free water) and 10 μL of viral RNA extracted from wastewater samples. All the quantitative assays were performed in sextuplicate. All the RT-qPCR reactions were performed using Hard-Shell 96-well PCR plates (Bio-Rad, USA) sealed with microseal “B” PCR Plate Sealing Film (Bio-Rad, USA). A non-template control (NTC) was included on each plate by replacing RNA with 10 μL of RNase-free water (Thermo Fisher Scientific, USA). Positive controls provided in the kit were also included for quality control.

Genomic library preparation and sequencing

Samples from the WWTP and from the municipality of A Coruña with lower Cq values were selected weekly for genomic library construction. cDNA synthesis was prepared from 11 μL of RNA using the SuperScript IV Reverse Transcriptase (Thermo Fisher Scientific, USA) and then amplified with ARTIC primer set v3 (Artic-network 2020). The resulting amplicons were mixed and purified with AMPure XP magnetic beads (Beckman Coulter, USA), and concentration was quantified with Qubit dsDNA HS Assay kit (Invitrogen, USA), attending manufacturer’s instructions. A total of 100 ng of DNA were taken into library preparation using the Illumina DNA Prep kit (Illumina, USA). Library concentration was measured using the Qubit dsDNA HS assay kit (Invitrogen, USA), and library validation and mean fragment size was determined using a Bioanalyzer DNA Analysis Kit (Agilent Technologies, USA). After dilution to 4 nM, libraries were pooled, denatured, diluted to 10 pM, and finally sequenced at read length 2 × 150 bp using the MiSeq Reagent Kit v2 (300 cycles) (Illumina, USA) and a MiSeq platform (Illumina, USA).

Measures of variables

Viral load was determined using the human 2019-nCoV RNA standard from European Virus Archive Global (EVAg). Calibration standards were prepared by diluting a stock solution of EVAg (10,000 copies/μL) in RNase-free water (Thermo Fisher Scientific, USA) within the range of 5–500 copies/μL for SARS-CoV-2. For each assay, a linear regression fit was established between log10 SARS-CoV-2 copy number and the Cq (quantification cycle) values for the N gene. A linear fit (y = mx + b) was obtained where y is SARS-CoV-2 RNA copies per L, x is the RT-qPCR Cq value, m is the slope, and b is the y-intercept. The calibration was included in every RT-qPCR run. Viral load in wastewater samples was calculated by using the slope and the y-intercept of the corresponding standard curves and the Cq values from RT-qPCR data.

The limit of detection (LoD) and the limit of quantification (LoQ) were determined from 20 replicates of the standard curves randomly chosen. LoD was defined as the lowest RNA concentration where all replicates were detected. LoQ was defined as the lowest concentration where the relative standard deviation (RSD) was < 30%. The standard curve parameters for the RT-qPCR assay such as the y-intercept (b), the slope (m), the linearity (R2), and the amplification efficiency (E) were also calculated following the recommendations by the MIQE guidelines (Bustin et al. 2009).

Models and data analysis procedure

Statistical smoothing methods

The use of nonparametric regression models for data smoothing is necessary to describe and understand the evolution of the pandemic. The fitting of these models is useful to identify outbreaks days from observing variables such as the viral load. The nonparametric approaches applied in this work were the Generalized Additive Models (GAM) (Wood 2017), the Locally Estimated Scatterplot Smoothing (LOESS) (Cleveland 1979), kernel regression (Wand and Jones 1994), local polynomial regression (Fan and Gijbels 1996 ), and the local bandwidth nonparametric regression (Herrmann and Maechler 2021).

GAM models are used given their flexibility to reproduce complex trajectories without involving experimental error fitting as a function of linear and smooth effects of the explanatory variables on the response (Vallejo et al. 2022). Thus, the expression Y = β0 + β1 X + s(T) + ε shows that the response variable \(Y\) can be expressed as a function of a linear effect of \(X\) and a smooth effect on the predictor \(T\), with \({\beta }_{0}\) and \({\beta }_{1}\) being the parameters defining the linear effect, \(s\left(T\right)\) the smooth effect, and \(\varepsilon\) the random error. In the present case, the smooth effect of time on viral load and the number of identified infected persons have been estimated, \(E(\mathrm{viral load}|\mathrm{tim}e)=s\left(\mathrm{time}\right)\) and \(E({N}^{\underset{\_}{^\circ }} \mathrm{of infected persons}|\mathrm{time})=s\left(\mathrm{time}\right)\) by the fitting of spline bases whose expression can be summarized by \(s\left(t\right)=\sum_{k=1}^{L+\mathrm{degree}}{\beta }_{k}{\phi }_{k}\left(t\right)={\beta }^{^{\prime}}\phi\), defined by a weighted sum of L (number of knots in the time interval) + degree (degree of the splines) elements, in which \({\phi }_{k}\left(t\right)\) are the spline functions. In this work, thin plate regression splines (fewer parameters to optimize), cubic regression splines with shrinkage (that include a penalty on the second derivative), and Gaussian process smooths (high flexibility) are used (Wood 2017).

A popular fitting alternative was the use of the local LOESS regression method (Cleveland 1979) incorporating two regression or local smoothing models, i.e., kernel and linear local polynomial. Kernel regression is also applied. Assume that the regression function is\(r\), and \({Y}_{i}\)(with \(i = 1, ..., n\)) are the \(n\) observations of the \(Y\) response variable. Thus, \(Y\) can be expressed as a function of the independent variable: \({Y}_{i}= r({t}_{i}) + {\varepsilon }_{i}\), with \({\varepsilon }_{i}\) the independent mean-zero error terms and \({t}_{i}\) the values of the \(t\) design variable where the model is evaluated. In this case, the regression function can be estimated using kernel methods, such as Gasser and Müller method (Gasser et al. 1991),\(\widehat{r}(t;h) = \frac{1}{h} \sum_{i=1}^{n}\left({g}_{i}\left(t;h\right)\cdot {Y}_{i}\right)\), where\({g}_{i}\left(t;h\right)=\underset{{s}_{i}-1}{\overset{{s}_{i}}{\int }}W\left[(t-u)/h\right]du\), si = (ti + ti+1)/2 with\({s}_{0}=0\),\({s}_{n}=1\), \(W\) is the kernel function (Gaussian, Epanechnikov, etc.), and h is the smoothing parameter, also called bandwidth. The latter is the expression used in the lokern package (Herrmann and Maechler 2021; R Core Team 2022), which also includes a plug-in method to estimate locally the bandwidth (Herrmann 1997). The optimal bandwidths are obtained by minimizing an estimate of the asymptotically optimal mean squared error (Herrmann and Maechler 2021).

On the other hand, the local linear regression estimator is also applied (Francisco-Fernández et al. 2015). Consider \({\{({t}_{i},{Y}_{i})\}}_{i=1}^{n}\) the observed values of a curve of viral load or number of infected persons, \({Y}_{i}\), evaluated at different times, \({t}_{i}\), with \(i=\mathrm{1,2},\dots ,n\), thus the he response variable can be expressed as a sum of the regression function \(r({t}_{i})\) as follows, \({Y}_{i}=r({t}_{i})+{\varepsilon }_{i}\), with \(1\le i\le n\) and \({\varepsilon }_{i}\) are the random errors. Of course, the aim is smoothing the data to reduce the noise of the experiment. This is done by fitting the data via local polynomial regression with a \(p\)-degree polynomial as shown by the expression \({\widehat{r}}_{h}(t)={{e}^{^{\prime}}}_{1}{({{X}^{^{\prime}}}_{t}{U}_{t}{X}_{t})}^{-1}{{X}^{^{\prime}}}_{t}{U}_{t}Y\), with \({e}_{1}={(\mathrm{1,0},\dots ,0)}^{^{\prime}}\), \(Y={({Y}_{1},\dots ,{Y}_{n})}^{^{\prime}}\), \({X}_{t}=\left[\begin{array}{c}1\\ \vdots \\ 1\end{array}\begin{array}{c}\left({t}_{1}-t\right)\\ \vdots \\ \left({t}_{n}-t\right) \end{array}\begin{array}{c}\cdots \\ \\ \cdots \end{array}\begin{array}{c}{\left({t}_{1}-t\right)}^{p}\\ \vdots \\ {\left({t}_{n}-t\right)}^{p}\end{array}\right]\) and \({U}_{t}=diag\{W(({t}_{1}-t)/h),\dots ,W(({t}_{n}-t)/h)\}\), while \(W(\cdot )\) is a kernel function, \(p\) the degree of the local polynomial (\(p=1\) for the local linear), and \(h\) the bandwidth. The bandwidth \(h\) is crucial to obtain a proper fit. In the specific case of the viral load fitting by a linear local polynomial, the bandwidth is locally estimated by using cross-validation (Vieu 1991).

The three alternatives for smoothing have provided a way to properly estimate the trend of viral load of COVID-19 in wastewaters. For the application of these techniques, the following R packages (R Core Team 2022), mgcv (Wood 2017), lokern (Herrmann and Maechler 2021), sm (Bowman and Azzalini 2021), and kernsmooth (Wand 2021), among others, have been used.

Estimation of the real number of infected people from SARS-CoV-2 viral load in wastewater

Statistical models previously developed by our team and described in detail in Vallejo et al. (2022) were used to estimate the positive cases in the metropolitan area of A Coruña served by the WWTP using the proportion of cumulative positive cases in the same area from the viral load data detected in wastewater. The population reported for each sampler has been obtained from the INE (Spanish National Statistics Institute). More concretely, data by census sections was first obtained from the statistics of the Continuous Register Statistics on January 1, 2019 (INE 2019), the last data available, and then, the census sections’ data corresponding to the sewage subnetwork of the sampler was aggregated.

Estimation of the frequency of SARS-CoV-2 variants

To estimate the frequency of the SARS-CoV-2 VOCs, we have previously implemented a maximum likelihood model that is described in detail in López de Ullibarri et al. (2023).

Bioinformatic analysis

Single-nucleotide variants (SNVs) and indels were identified with iVar (Grubaugh et al. 2019). We aligned the paired-end reads to the reference sequence MN908947.3 from Wuhan using the BWA-mem tool (Li 2013) and sorted them using SAMtools (Danecek et al. 2021). Then, iVar trim was used to soft-clip the primer sequences and low-quality bases and remove reads shorter than 30 bp. Then, iVAR consensus was applied, with a minimum frequency threshold of 0.5. In addition, iVar functions getmasked and removereads were used to remove the reads corresponding to the indexes of the mismatched primers. SAMtools mpileup and iVar variants were used to identify single nucleotide changes and indels. A minimum frequency threshold of 0.0001 and a minimum depth of 1 were set for running iVar. The minimum base quality score threshold was left as default, 20. Then, iVar’s output was processed in R software (R Core Team 2022) to detect typical mutations of specific SARS-CoV-2 variants, which were obtained from the website http://outbreak.info.

Open access data sharing

A public website was developed to show the current epidemiological situation over time in A Coruña (Spain), available at https://edarbens.es/covid19/. The Leaflet library (Agafonkin 2020) was used to include maps showing the historical viral load (copies of viral RNA per L) in the different districts over time. These maps use viral load data included as JSON on the website. Daily estimations of the viral load are also publicly available on the website. COVIDBENS is included in the NORMAN SCORE “SARS-CoV-2 in sewage” (SC2S) database (NORMAN 2021), an international platform for rapid, open access data sharing about wastewater pathogen surveillance.

Results

RT-qPCR standard curve characteristics

LoD and LoQ values of the RT-qPCR assay were determined using 20 replicates of serial dilutions of EVAg standard ranging from 5 to 500 RNA copies/µL. Cq values and detection rates are shown in Table S1 (Online Resource 1), showing that 10 copies/µL is the LoD and 5 copies/µL is de LoQ. The slope of the standard curve for the N gene was − 3.1709, the y-intercept was 41.134, the R2 value was 0.9918, and the amplification efficiency was above 106.71% (Online Resource 2).

Evolution of SARS-CoV-2 viral load over time

For the present study, viral load was analysed over 22 months in the metropolitan area of A Coruña (Spain). The COVIDBENS monitoring program collected a total of 863 samples of sewage from the Bens WWTP and from the municipalities of A Coruña, Arteixo, Culleredo, Oleiros, and Cambre. RT-qPCR results of viral load (copies of viral RNA per L) and Cq values obtained for the N gene are given in Online Resource 3.

SARS-CoV-2 RNA was detected in 96.5% of the wastewater samples. Detection rates of the viral RNA of the different sampling points are represented in Fig. 2. Viral load calculated from the N gene ranged from 1.30 × 103 to 5.87 × 106 copies/L (Table 1). The detection of SARS-CoV-2 in wastewater in the metropolitan area of A Coruña, including the five municipalities, collected in the WWTP of Bens (Fig. 3) showed six viral load waves in wastewater after the first COVID-19 outbreak reported in March 2020 (Vallejo et al. 2022). After the lockdown and the imposition of new restriction measures by the government the lowest viral load in wastewater was obtained, failing to amplify between June and July 2020. Nevertheless, in July 2020 the viral load increased again beginning the first wave, obtaining a maximum of 218.755 copies per L. The second wave covered the period from September to November 2020, where a maximum of 831.309 copies per L was obtained. After this second wave, viral load in wastewater increased again making way for the third wave that began with the emergence of the Alpha variant (B.1.1.7) in January 2021, reaching a maximum of 1.276,329 copies per L. Then, the viral load decreased again until the emergence of the fourth wave, where a peak between March and April 2021 was detected. Although the viral load values remained high, a fifth wave was observed from May to September 2021, where viral load reached around 3 million copies per L due to the appearance of the Delta variant (B.1.617.2). The last wave reported in this study is represented by the progressive increase in viral load from November 2021 to February 2022 due to the irruption of Omicron variant (B.1.1.529), obtaining a maximum of 5.102.970 copies per L. The evolution of the viral load measured separately in the different municipalities, which started some months later than in the whole area, followed a similar trend (Fig. 4), being consistent with that observed at the WWTP of Bens.

Fig. 2
figure 2

Detection rates of SARS-CoV-2 in wastewater samples of the different sampling points located in the metropolitan area of A Coruña (Spain). N represented the number of collected samples

Full size image
Table 1 Summary of SARS-CoV-2 monitoring results at the different sampling points in the metropolitan area of A Coruña (Spain) from June 2020 to March 2022
Full size table
Fig. 3
figure 3

Evolution of the SARS-CoV-2 viral load measured at the Bens WWTP (A Coruña, Spain). The transition to each new variant is indicated with grey shadows

Full size image
Fig. 4
figure 4

Evolution of the SARS-CoV-2 viral load in samples obtained from sewage from the municipalities of A Coruña, Arteixo, Cambre, Culleredo, and Oleiros from December 19, 2020 to March 17, 2022. The transition periods to new variants are indicated with grey shadows

Full size image

COVIDBENS as an early warning system

SARS-CoV-2 viral load measured at the Bens WWTP was plotted and compared with the number of active clinical cases reported by the health authorities in A Coruña (Fig. 5). The comparison revealed that an increase in viral load in wastewater clearly preceded a subsequent increase in clinical cases in each of the successive pandemic waves. Official COVID-19 cases reported by the health authorities are given in Online Resource 4. Remarkably, increases in viral load corresponding to the different epidemic waves were detected in wastewater before the clinical outbreaks were reported (Fig. 5) with 8–36 days of anticipation (Table 2).

Fig. 5
figure 5

Representation of the early warning system as assessed by COVIDBENS. Smoothed viral load measured at the Bens WWTP (red line), official number of active cases (blue line), reported by the health authorities (XUGA, Spain), and percentage of vaccinated population (green line) in the metropolitan area of A Coruña (Spain). Vertical dotted lines highlight increases of viral load (red) and clinical cases (blue). The transition periods to new emerging variants are indicated with grey shadow

Full size image
Table 2 SARS-CoV-2 detection results and prediction of COVID-19 clinical cases (lead time) in this study (A Coruña, Spain) from June 2020 to March 2022, and in other studies
Full size table

The monitorization of SARS-CoV-2 mutations in wastewater was carried out from November 2020 to March 2022. Moreover, the evolution of variants’ frequency in the local population over time was determined from May 2021 to March 2022 (Fig. 6) using the statistical method previously described by our team (López de Ullibarri et al. 2023). We detected the Alpha variant in the wastewater samples on December 16, 2020, 42 days before it was detected in the clinical samples (January 27, 2021) (Fig. 6), being the longest lead time achieved here (Table 3). Thus, Alpha was the predominant variant during the third viral load wave. Similarly, we detected the Delta variant in wastewater on May 18, 2021, 30 days before the appearance of clinical cases (June 17, 2021). Delta became predominant during the fifth viral load wave. After that, the Delta abundance in wastewater decreased progressively until it was replaced by the Omicron variant in January 2022. Finally, we could also anticipate the emergence of the Omicron lineage BA.2, detected in wastewater on January 18, 2022, 27 days before its identification in clinical samples on February 14, 2022.

Fig. 6
figure 6

Evolution of the SARS-CoV-2 variant frequency over time based on mutations detected in wastewater samples from the metropolitan area of A Coruña (Spain) from May 2021 to March 2022 using statistical models previously described. Arrows indicate the first clinical case reported by the health system infected with each variant (data recovered by the Galician Health Service, SERGAS)

Full size image
Table 3 SARS-CoV-2 variants predicting results in this work (A Coruña, Spain) from November 2020 to March 2022, and in other studies
Full size table

Evolution of the estimated number of infected people over time

Statistical models previously developed by our team (Vallejo et al. 2022) were leveraged to estimate the number of people infected by SARS-CoV-2 over time in the successive epidemic waves (Online Resource 5), including symptomatic and asymptomatic people. The SARS-CoV-2 viral load in wastewater (Fig. 3) and the estimated number of infected people followed similar trends (Fig. 7) and peaks coincided with the emergence of SARS-CoV-2 variants.

Fig. 7
figure 7

Evolution of the number of people infected by SARS-CoV-2 in the metropolitan area of A Coruña (Spain) during the pandemic from July 2020 to March 2022. Data have been estimated using COVIDBENS statistical models (Vallejo et al. 2022)

Full size image

Discussion

Although many efforts have been made by researchers all over the world to detect asymptomatic people infected with SARS-CoV-2, using nasopharyngeal PCR or antigen and antibodies tests, to date, we are still not able to identify the entire population infected and, therefore the real magnitude of the epidemic is difficult to understand. Wastewater monitoring stands as a powerful tool to monitor the COVID-19 epidemic in the entire infected population, including symptomatic and asymptomatic people.

SARS-CoV-2 was absent in wastewater during the period June-July 2020, which can be explained by the efficacy of the hard restriction measures implemented from mid-March to the end of June. But, unfortunately, at the beginning of July the first new epidemic outbreak after restrictions occurred, followed by a long succession of epidemic waves. During the period of the present study, the maximum viral load detected in wastewater samples ranged from 105 to 106 RNA copies per L, similar to that reported in other studies (Gerrity et al. 2021; Padilla-Reyes et al. 2022; Pillay et al. 2021; Saththasivam et al. 2021). However, Yaniv et al. (2021) showed higher viral load values, around 107 RNA copies per L, which could be related to the length of their study (4 days) and to the type of sample collecting (8 h composite samples). On the contrary, Kumar et al. (2021) reported a low viral load range, 104 RNA copies per L, probably due to the use of grab samples, which, as reported by Rafiee et al. (2021), makes difficult to follow the evolution of viral load or to perform statistical approaches. Wastewater epidemiology is influenced by many factors such as the SARS-CoV-2 shedding pattern of each variant, the population size of each location, the physico-chemical characteristics of the wastewater, and others such as the commercial activity, or the international trade within and between countries, which induces variability in viral load data. Accordingly, the shedding of SARS-CoV-2 RNA in the human body depends on the stage of the disease in which the patient is. Wu et al. (2022) demonstrated an early SARS-CoV-2 shedding peak before symptom onset with a subsequent progressive reduction in RNA load in the last stages of the disease. The age also affects viral shedding, in such a way that a higher viral load is detected in older people (Vellas et al. 2020). Moreover, it can be established that the recovery of COVID-19 patients implies less shedding. Also, the total population covered by the sampling area is a crucial factor in WBE since the larger the population, the greater the flow of wastewater, thus diluting the viral RNA that reaches the WWTP (Wilder et al. 2021). It also should be taken into account that tourism and travel, which significantly change during summer, vary the population size. In addition, wastewater properties such as pH, temperature, humidity, ammonia or solid composition, or chlorination are directly associated with SARS-CoV-2 viral load measured in sewage samples (Ahmed et al. 2020; Amoah et al. 2022; Zhang et al. 2020). SARS-CoV-2 detection differs between the solid and the liquid phase of the wastewater influent, reporting a higher viral load in the liquid fraction (Weidhaas et al. 2021). In our previous work (Vallejo et al. 2022), we suggested that a high flow rate due to the rainfall and long transit times through the sanitary network until the WWTP, reduced SARS-CoV-2 concentration in wastewater, because the viral RNA was exponentially degraded during transportation. Therefore, the population living near the WWTP is probably the main contributor to the viral load detected in wastewater. The dilution effect of rainfall reduces SARS-CoV-2 RNA levels in wastewater. Saingam et al. (2023) used chemical indicators (phosphate and ammonia) to calibrate this dilution effect improving the correlation between viral load in wastewater and clinical data over time. Lastly, one of the factors that influence the amount of SARS-CoV-2 in wastewater and the detection of new mutations is international trade. This factor is a good indicator of the spread of the COVID-19 disease, since the most active areas in trade have greater population mobility and more group activities, increasing interpersonal contacts. Thus, social interactions are the main source of the spread of SARS-CoV-2 in society and it is correlated with confirmed clinical cases (Bontempi and Coccia 2021). In addition, international trade promotes the contact with foreign populations and consequently increases the risk of SARS-CoV-2 transmission between communities (Bontempi et al. 2021). These indicators help in the design of new effective policy responses to reduce the impact of future pandemics on society.

COVIDBENS was able to anticipate up to 36 days’ increases in clinical cases, or epidemic waves. Similar WBE studies resulted in 7–24 days of anticipation (Claro et al. 2021; Kumar et al. 2021; Robotto et al. 2022; Sangsanont et al. 2022). Remarkably, in Barcelona, Chavarria-Miro et al. (2021) reported an outbreak 41 days before the first clinical sample was detected. They performed different RT-qPCR assays targeting 5 SARS-CoV-2 genes instead of only one, increasing the robustness of the data. In other studies, despite using 24-h composite samples over a long surveillance period, no lead time was observed, as in the study reported by Reynolds et al. (2022). The ability to anticipate to early outbreaks depends on the sampling strategy used, the frequency of sampling and on the smoothing statistical methods used, among other factors. It has to be taken into account that, in general, limited data or punctual analyses do not allow establishing a reliable relationship between viral load and clinical cases using statistical models. COVIDBENS recruited samples at the Bens WWTP including all the municipalities of the metropolitan area of A Coruña, but also monitored the wastewater from specific locations, which was essential to detect local outbreaks. Regarding the sampling frequency, we collected samples twice or three times per week, during the period of the present study, depending on the evolution of the epidemic. After a preliminary study, we chose 24-h composite sampling as the best scheme for monitoring the epidemic in the metropolitan area of A Coruña.

A large proportion of people infected by SARS-CoV-2 have no symptoms, so they are not usually detected and thus reported by the health system. The statistical model used in this study (Vallejo et al. 2022) allowed us to estimate the proportion of people infected in the metropolitan area of A Coruña. However, this model was developed when lineage B.1.177 was prevalent in Spain. After that, several new VOCs have emerged with potentially larger viral shedding and consequently higher load in the wastewater. This fact requires an adaptation of the statistical models to the new situation of higher viral load excreted per person as new variants appear. Several studies confirmed the longer duration of Alpha variant in the respiratory tract of infected people resulting in higher RNA load compared to previous lineages (Calistri et al. 2021; Lyngse et al. 2021). Similarly, both Delta and Omicron variants seem to result in higher RNA load per patient than earlier variants (Riediker et al. 2022). Other authors suggested that the Omicron was more transmissible and contagious than Delta, but nevertheless they conclude that this variant does not generate a higher viral load per person than Delta (Chen et al. 2022; Migueres et al. 2022; Sentis et al. 2022). This explains the higher viral load values observed in wastewater in January 2022 due to the rapid spread of Omicron in the population. Therefore, epidemiological surveillance models for SARS-CoV-2 in wastewater must consider these differences in the levels of excretion of VOCs. Results of the estimated number of people infected by SARS-CoV-2 until January 2021 are 90% reliable as described by Vallejo et al. (2022), but from then on, a specific model adjustment would be necessary for each new variant emerged over time.

A relevant problem in the context of the COVID-19 pandemic was the monitorization of the progression of SARS-COV-2 variants. The most used classification tools, such as pangolin, were designed for clinical samples in which only one lineage predominates, but in case of wastewater samples, several variants coexist at the same time, which makes the bioinformatics analysis harder. SARS-CoV-2 genome sequencing from wastewater allows the detection of emerging variants far before they appear in clinical samples. We reported the Alpha, Delta, and Omicron variants up to 42 days earlier than the clinical cases detected by the regional health system. This anticipation demonstrates the great potential of genomic surveillance in wastewater samples. Ideally, for best WBE, variant detection and viral load monitoring in wastewater over time should be investigated in a complementary way (Galani et al. 2022; Johnson et al. 2022; Masachessi et al. 2022; Rubio-Acero et al. 2021). Few studies reported similar lead times (Joshi et al. 2022; Vo et al. 2022) while other achieved shorter lead times, around 7 days (Kirby et al. 2022). The genomic surveillance method used by Kirby et al. (2022) does not allow knowing if all the variant-associated mutations are present in a single genome, making difficult to confirm the presence of a specific variant in the wastewater sample. However, the statistical methods previously developed by our team (López de Ullibarri et al. 2023), used in the present study, allowed the determination of the percentage of each variant in the population over time from mutations data obtained from pool samples such as wastewater. More and more studies are implementing bioinformatic tools and statistical approaches for the early detection of low-frequency variants and the quantitative monitoring (Jahn et al. 2022). Karthikeyan et al. (2022) developed software for SARS-CoV-2 variants differentiation from wastewater samples. In addition, it was shown that the combined use of arctic v4 primers and the Maxwell® RSC Enviro Wastewater TNA Kit (TNA) increases the sensitivity and genome coverage of the virus sequencing (Girón-Guzmán et al. 2023).

The irruption of the Delta variant in the summer of 2021 had a much greater impact on viral load in wastewater than the Alpha variant had, probably due to its mutations in the spike protein which highly impact in its transmissibility and immune evasion (Tian et al. 2021). However, the number of COVID-19 reported cases during the Delta wave was lower, reflecting the vaccination effectiveness. The transition to the Omicron variant corresponds to the last viral load wave reported in this study. Despite the high vaccination rate in A Coruña (above 90%), the higher transmission rate of Omicron compared to Delta (Allen et al. 2022; Kumar et al. 2022) and its immune escape ability (Zhang et al. 2021) probably explain a larger viral load in the wastewater, in contrast to the low rate of hospitalized people (data not shown), which indicates the benefit of vaccines. Omicron continues to accumulate mutations, emerging several lineages which have rapidly spread globally during 2022 (Parums 2022). The BA.2 sublineage became the dominant variant at the end of this monitoring program due to its higher transmissibility and its capacity to reduce the effect of vaccination (Lyngse et al. 2022). Currently, the XBB.1.5 Omicron, detected for the first time in the USA, became predominant in Spain. The great impact of the VOCs in the population and their capability of reducing the protective effect of the vaccines increased the transmissibility of the infection, reporting a higher proportion of COVID-19 cases. However, hospitalization and mortality in COVID-19 cases were significantly reduced as the percentage of vaccinated people increased. This demonstrates the effectiveness of mass vaccination campaigns to contain the COVID-19 pandemic.

SARS-CoV-2 surveillance in wastewater in A Coruña helped in decision-making at both the political and hospital levels. Both the regional and local government used COVIDBENS data to implement the corresponding restriction measures, and hospitals were able to anticipate and prepare for an increase of cases. In addition, important companies have followed our data to plan their production according to the COVID-19 pandemic, suggesting that industry, among other sectors, can benefit from the implementation of these epidemiological models. Nevertheless, in addition to the implementation of an effective WBE system, it is essential to have effective and preventive strategies against future pandemic threats, not only focused on the health level but also on the social, environmental, and institutional levels. During the COVID-19 pandemic, it has been shown that measures such as the use of masks or lockdown considerably stopped the transmission of SARS-CoV-2 in the population. Other necessary preventive measures are the control of the transport and trade of pathogenic microorganisms, the reduction of air pollution levels in big cities, and transparent, responsive, and effective governance (Coccia 2021). Strategies such as contact tracing and clinical testing to diagnose SARS-CoV-2 infections in the population would reduce the transmission of the disease, minimizing the number of infected people and avoiding the greatest possible number of deaths, especially in the initial stages of a pandemic situation when effective drugs or vaccines are not yet available (Benati and Coccia 2022). In addition, the implementation of a rapid vaccination plan at the beginning of the pandemic, the negative impact of the pathogen on society would be reduced and future outbreaks and similar pandemics would be prevented (Coccia 2022a). For this reason, countries have to increase the R&D investment for the development of new effective vaccines and drugs. Coccia (2022b) developed the r (resilience) and p (preparedness) indexes, which measure the ability of countries to face the COVID-19 crisis based on their level of support for vaccine development and their mortality reduction policies. It is an effective method that could be implemented in prevention strategies for future similar pandemics. Even so, new improvements and scientific advances are still needed on WBE, not only focused on SARS-CoV-2 but also on other pathogenic microorganisms with pandemic potential (Núñez-Delgado et al. 2021).

Even though COVIDBENS monitoring program ended in March 2022, SARS-CoV-2 wastewater, we strongly believe that wastewater surveillance should continue for monitoring the COVID-19 pandemic. WBE is an effective tool with various applications that support public health actions (Prado et al. 2023): (i) provides information about the transmission of SARS-CoV-2 and other potentially dangerous pathogens; (ii) detects early outbreaks; (iii) allows monitoring the effectiveness of the measures imposed by governments and local authorities to contain the viral transmission; (iv) alerts about the emergence of new pathogens, such as Crimean Congo Haemorrhagic Fever virus or Monkeypox virus, recently considered in Europe; (v) serves to monitor antimicrobials resistance in the community; and (vi) helps political authorities to be prepared for future pandemics. For these reasons, the European Commission recommended the inclusion of WBE in the national detection strategies (EC 2021b).

Limitations of the present study

Our study has some limitations that must be taken into account. First, we did not consider the dilution rate in wastewater. Flow normalization is commonly used when reliable flow data are available, but in the absence of these, alternative methods such as electrical conductivity and crAssphage can be very effective (Langeveld et al. 2023). Secondly, we assumed no variation in the shedding pattern along the surveillance period, but the different VOC excretion levels in feces affect both the early warning condition of the WBE for COVID-19 and the estimated infected people results. This could be related with a decrease in the ability to anticipate the emergence of Omicron with respect to Delta. As we described before, our statistical model to estimate the people infected with SARS-CoV-2 was only designed for the original variant (B.1.177), so it must be adjusted for new each new emerging variant over time. Thirdly, WBE has a low sensitivity and recovery of the SARS-CoV-2 detection in sewage due to the complexity of the wastewater samples, so new methods are needed to improve the sensitivity and virus recovery such as the EM-VIP-Mag-RT-qPCR method (Kumblathan et al. 2023). In addition, due to the distribution of the sewerage network, the specific location of the sewers installed in the five municipalities of the metropolitan area of A Coruña do not cover the entire population, which may affect the WBE system efficiency, given that we lost a small part of the infected population. A deeper redesign would be necessary to resolve this issue. Lastly, the use of automatic samplers is highly expensive, so it is not valid for all types of regions. Studies have demonstrated the efficacy of passive sampling in WBE in small areas using Moore’s swabs, which are more sensitive, precise, and cost-effective, so they could be implemented in low-resource areas where clinical testing is scarce (Cha et al. 2023).

Conclusions

This work demonstrates that SARS-CoV-2 viral load monitoring in wastewater is a very efficient way to report on COVID-19 evolution in the community, anticipating new outbreaks long before the health system, and therefore allowing the public authorities to take reliable decisions, saving precious time and public resources. SARS-CoV-2 monitoring in wastewater detects increases in viral load before clinical cases appear being able to predict epidemic waves in the population. Viral load data combined with statistical models allows estimating the total number of infected people (symptomatic and asymptomatic), improving significantly the WBE strategy. In addition, sequencing of SARS-CoV-2 in wastewater allows a real-time epidemiological surveillance of SARS-CoV-2 mutations and combined with statistical models, can anticipate the appearance of emerging VOCs before clinical testing. COVIDBENS has served as an effective early warning system for anticipating new outbreaks and variants of concern in the metropolitan area of A Coruña and, it also provided a public health service open to all citizens with an important social impact throughout the entire pandemic. Nations should implement the WBE system in their prevention strategies for future pandemic crises, which also include social, environmental, and health policies, to reduce the socioeconomic impact of the virus or other potentially dangerous pathogens.