Dynamics of SARS-CoV-2 Alpha (B.1.1.7) variant spread: The wastewater surveillance approach.

Wastewater based epidemiology (WBE) offers an overview of the SARS-CoV-2 variants circulating among the population thereby serving as a proper surveillance method. The variant of concern (VOC) Alpha was first identified in September 2020 in the United Kingdom, and rapidly became dominant across Europe. Our objective was to elucidate the Alpha VOC outcompetition rate and identify mutations in the spike glycoprotein (S) gene, indicative of the circulation of the Alpha VOC and/or other variants in the population through wastewater analysis. In the period covered by this study (November 2020-April 2021), forteen wastewater treatment plants (WWTPs) were weekly sampled. The total number of SARS-CoV-2 genome copies per L (GC/L) was determined with a Real-Time qPCR, targeting the N gene. Surveillance of the Alpha VOC circulation was ascertained using a duplex RT-qPCR, targeting and discriminating the S gene. Our results showed that in a period of 6 weeks the Alpha VOC was present in all the studied WWTPs, and became dominant in 11 weeks on average. The outcompetition rates of the Alpha VOC were estimated, and their relationship with different parameters statistically analyzed. The rapid spread of the Alpha VOC was influenced by its initial input and by the previous circulation of SARS-COV-2 in the population. This latter point could be explained by its higher transmissibility, particularly advantadgeous when a certain degree of herd immunity exists. Moreover, the presence of signature mutations of SARS-COV-2 variants were established by deep-sequencing of the complete S gene. The circulation of the Alpha VOC in the area under study was confirmed, and additionally two combinations of mutations in the S glycoprotein (T73A and D253N, and S477N and A522S) that could affect antibody binding were identified.

The corresponding author confirm that all authors agree with the author lineup, the order of the names, and the content of the manuscript.

Introduction

In September 2020, the SARS-CoV-2 variant of concern (VOC) Alpha (B.1.1.7) emerged in Southeast England, and rapidly outcompeted other variants and became dominant across the UK and Europe (https://assets.publishing.service.gov.uk/government/uploads/system/uploads/attachment_data/file/959426/Variant_of_Concern_VOC_202012_01_Technical_Briefing_5.pdf; (Davies et al., 2021). Although spreading of the different SARS-CoV-2 lineages relies on their relative fitness, founder effects, as well as on the social behavior of the infected patients (LauringHodcroft, 2021), the rapid propagation of the Alpha variant has been mostly related with its higher transmissibility, i.e. fitness, compared with the previously dominant variant B.1.177 (Davies et al., 2021; Frampton et al., 2021; Volz et al., 2021).

Sequencing of clinical specimens revealed that introduction of the Alpha variant in Spain occurred through several independent establishments, starting late December 2020 (Buenestado-Serrano et al., 2021; Perramon et al., 2021), and thereafter progressively spread across the country. However, sequencing of clinical samples covers only 1–2% of cases (https://www.mscbs.gob.es/profesionales/saludPublica/ccayes/alertasActual/nCov/documentos/Integracion_de_la_secuenciacion_genomica-en_la_vigilancia_del_SARS-CoV-2.pdf), mostly from big cities, weakening the analysis of the spatiotemporal evolution of VOCs.

Surveillance of SARS-CoV-2 in sewage may be very useful for the study of the dynamics of the spread of VOCs in a given geographic region, since it enables detection of genomes from symptomatic and asymptomatic patients and also permits the analysis of many wastewater treatment plants (WWTP) including highly and low populated areas (Heijnen et al., 2021; Lee et al., 2021). In consequence, wastewater analysis provides information on the first introduction events of the variants, their spatiotemporal distribution and the overall competition between them.

In this study, we present the evolution of the Alpha VOC in the Autonomous Community of Catalonia in Spain, through the analysis of 14 WWTPs located across the Community.

Materials and methods

Wastewater sampling

Twenty-four hour composite 1-l samples of influent wastewater from 14 WWTPs (Table 1 ), covering a population of 2,455,668 inhabitants (32% of the total Catalan population) in an area of 1,757 Km2, were weekly collected within the framework of the Catalan Surveillance Network of SARS-CoV-2 in Sewage (https://sarsaigua.icra.cat/; http://doi.org/10.5281/zenodo.4147073; {Corominas et al., 2021 #1}). Samples were transported at 0°-4°C in a portable icebox, concentrated upon arrival, and analyzed the next day. The present study covers the period November 2020 to April 2021.

Table 1

Wastewater treatment plants (WWTPs) included in the study.

Wastewater Treatment Plant	Sample collection	Sewer Length (Km)/number of pumps	Population equivalents	Population covered	Population density (Inhabitants/Kmb)
Puigcerdà (PUI)	Compositea (24 h)	19/1	25,115	10,675	330
Berga (BER)	Compositea (24 h)	5/1	15,880	16,494	742
Riera de la Bisbal (VEN)	Compositeb (24 h)	38/17	60,779	50,062	837
Abrera (ABR)	Compositec (24 h)	30/7	69,376	65,160	949
Vilafranca del Penedès (VDP)	Composited (24 h)	33/18	67,763	67,381	1,300
Banyoles (BAY)	Compositea (24 h)	11/5	39,554	27,959	1,319
Martorell (MRT)	Compositec (24 h)	12/1	26,425	28,189	1,866
Rubí (RUB)	Composited (24 h)	13/0	102,148	168,470	2,159
Figueres (FIG)	Compositea (24 h)	15/5	79,139	52,048	2,262
Girona (GIR)	Compositea (24 h)	55/6	120,059	153,162	2,749
Granollers (GRA)	Compositea (24 h)	22/0	134,494	97,996	3,200
Montcada i Reixach (MIR)	Compositee (24 h)	30/0	248,261	269,784	3,926
Sabadell Riu Sec (SRS)	Compositef (24 h)	12/1	155,094	266,876	5,196
El Prat de Llobregat (PDL)	Compositea (24 h)	71/8	1,263,797	1,181,412	8,918

(Proportional flow every 20 min).

(Non-proportional flow every hour).

(Non-proportional flow every 20 min).

(Proportional flow every hour).

(Non-proportional flow every 15 min).

(Proportional flow every 2 h).

Wastewater concentration

Two hundred ml of wastewater were concentrated, using the aluminum hydroxide adsorption-precipitation method (Randazzo et al., 2020; WallisMelnick, 1967), to a final volume of 1–5 ml of phosphate buffered saline (PBS). As process control, 1.5 × 106 TCID50 units of the attenuated PUR46-MAD strain of the Transmissible Gastroenteritis Enteric Virus (TGEV) (Moreno et al., 2008) were seeded into each sample previously to the concentration step.

Nucleic acids extraction

Nucleic acids were extracted from 300 μl of concentrates using the Maxwell® RSC PureFood GMO and Authentication Kit (Promega), following the manufacturer's instructions. For each extraction run, a PBS negative control was included. To estimate virus recovery a previously described RTqPCR for TGEV quantification was applied (Vemulapalli and GulaniSantrich, 2009). Samples with a virus recovery ≥1% were considered acceptable. Our recoveries ranged between 1% and 98%, with an average of 24% ± 18%.

RTqPCR assays

Quantification of SARS-CoV-2 RNA in sewage samples was based on the N1 assay (US-CDC, 2020), which targets a fragment of the nucleocapsid gene. We employed the PrimeScript™ One Step RT-PCR Kit (Takara Bio, USA) and a CFX96 BioRad instrument.

A duplex TaqMan RTqPCR assay for allelic discrimination of the Spike gene (S), based on the detection of a deletion of six nucleotides (nt position 21,765–21,770; amino acid position 69–70) in the Alpha variant vs. the B.1.177 variant, was implemented (Carcereny et al., 2021). An 89-nt fragment of the S gene was amplified (primers For-S21708 5′-ATTCAACTCAGGACTTGTTCTTACCTT-3′ and Rev-S21796 5′-TAAATGGTAGGACAGGGTTATCAAAC-3′) and detected with two probes (69-70del Probe: 5′HEX-TTCCATGCTATCTCTGGGACCAATGGTACT BHQ1-3′ and 69-70in Probe: 5′FAM-TCCATGCTATACATGTCTCTGGGACCAATG BHQ1-3′) targeting the Alpha and B.1.177 amplicons, respectively, using the PrimeScript™ One Step RT-PCR Kit (Takara Bio, USA). Although the 69–70 deletion (69-70del) is not exclusive of the Alpha variant, in the region and period under study, other variants bearing this deletion were very uncommon (Carcereny et al., 2021); https://canalsalut.gencat.cat/web/.content/_A-Z/C/coronavirus-2019-ncov/material-divulgatiu/informe-vigilancia-variants-genomiques-sarscov-2.pdf).

To monitor the presence of inhibitors, RTqPCR assays for each target included two wells of undiluted and two wells of ten-fold diluted RNA. Additionally, each RTqPCR assay included four negative control wells (two of nuclease-free water and two negative extraction controls). Standard curves for genome copies quantitation were performed using two commercial synthetic SARS-CoV-2 RNA controls (Twist Control 2, MN908947.3; and Twist Control 14, EPI_ISL_710,528). The limit of quantification (LOQ) for each specific target was determined by running series of dilutions of each target with 4–10 replicates per dilution, and was established in 28.45 genome copies per reaction (GC/rx), 27.61 GC/rx and 24.76 GC/rx for the N1, the Alpha variant and the B.1.177 S targets, respectively (Carcereny et al., 2021).

Cq values ≤ 40 were converted into genome copies/L (GC/L) using the standard curves and taking into consideration the volumes tested. Inhibition was established when differences in viral titers obtained using the undiluted and the ten-fold diluted RNAs were higher than 0.5 log10 in which case titers were estimated using only the diluted values. Otherwise, the average of the four values was applied.

For comparison analysis between different WWTPs and weeks, SARS-CoV-2 genome titers (GC/L) based on the N1 target quantitation, were normalized vs. the ammonium content (data provided by the WWTPs; see Table S1) dividing the GC/L by the NH4 + concentration in mg/L (Been et al., 2014). Additionally, to get SARS-CoV-2 genome titers per 100,000 inhabitants we used the following formula:

Finally, the log10 GC/100,000 inhabitants was further normalized by the population density of each WWTP, which was estimated taking into consideration the population density (https://www.idescat.cat/) of each of the served municipalities, and their percentage of contribution to the total WWTP flow volume (http://aca.gencat.cat/ca/laigua/infraestructures/estacions-depuradores-daigua-residual/).

The percentage of SARS-CoV-2 genomes bearing the 69-70del in the S gene out of the total S gene copies in each sewage sample was estimated using the formula:

In samples with one of the concentrations < LOQ, the percentage was calculated using the LOQ value of the assay. Samples with both concentrations < LOQ were not included in the dynamics analysis.

The rate of outcompetition of the Alpha variant in each WWTP was assessed by calculating the slope of the regression line between the variable time (in weeks) and the variable percentage of Alpha variant. Regression was calculated using SigmaPlot v. 11.0 (Systat Software Inc, CA, USA).

The log10 GC/100,000 inhabitants of the Alpha variant was obtained by applying the percentage of the 69-70del to the N1 log10 GC/100,000 inhabitants, which were previously normalized by ammonium content and population density as described above in this section.

COVID-19 cases in the areas covered by the WWTPs

The number of accumulated COVID-19 cases associated to each WWTP was calculated by summing up the number of cases of each of the municipalities covered by each of the WWTPs. Municipalities associated to each WWTP were obtained from http://aca.gencat.cat/ca/laigua/infraestructures/estacions-depuradores-daigua-residual/and the number of cases were obtained from the following public servers: http://governobert.gencat.cat/ca/dades_obertes/dades-obertes-covid-19/and https://biocomsc.upc.edu/en/covid-19.

The number of cases were expressed per 100,000 inhabitants and were also normalized by population density.

Full-length S gene deep sequencing

Samples with Cq ≤ 38 at the 10-fold diluted wells, were chosen for sequencing with the ARTIC Network protocol (https://artic.network/ncov-2019) using selected v3 primers (Integrated DNA Technologies) for genome amplification, and KAPA HyperPrep Kit (Roche Applied Science) for library preparation. Briefly, libraries were prepared as follows: (i) cDNA synthesis using random hexamers, (ii) target-specific PCRs of the complete S gene (14 amplicons in two separate pools) encompassing 40 cycles of amplification, and DNA purification using the Kapa Pure Beads (Roche), (iii) quantitation using the Qubit® dsDNA HS Assay Kit (Thermo Fisher Scientific, Waltham, MA, USA) and normalization to 1.5 ng/μl with Tris-HCl 10 mM, (iv) end repair, A-tailing and adaptor ligation using the Illumina Kapa Single-Indexed Adapter Kit (Set A + B), and DNA purification using the Kapa Pure Beads, (v) 8-cycle PCR amplification of libraries, and DNA purification using the Kapa Pure Beads and (vi) Qubit quantitation and normalization to 4 nM with Tris-HCl 10 mM. Libraries were sequenced using MiSeq Reagent Kit 600v3 cartridges on the MiSeq platform (Illumina). On average 50,000–70,000 reads per amplicon were obtained.

The bioinformatics analysis included the next steps: (i) MiSeq R1 and R2 paired ends were used to reconstruct each amplicon using FLASH1 program, setting a minimum of 20 overlapping bases and a maximum of 10% mismatches, (ii) reads not meeting minimum requirements were discarded. Reads with more than 5% of bases below a Phred score of Q30 were filtered out, (iii) reads were demultiplexed by matching primers, allowing a maximum of three mismatches, (iv) primers were trimmed at both read ends, (v) identical reads were collapsed to haplotypes, (vi) the frequency of each haplotype corresponded to the read counts, (vii) a fasta file per pool/primer/strand combination was generated, and reverse haplotypes were reverse-complemented, (viii) raw forward and reverse haplotypes were multiple aligned with MUltiple Sequence Comparison by Log-Expectation (Muscle), (ix) haplotypes common to both strands at abundances ≥0.1% were maintained and those below that value or unique to one strand were discarded and (x) haplotypes common to both strands and with frequencies ≥0.1% are selected for subsequent computations (Andres et al., 2020). Clean sequences were aligned against the Wuhan-Hu-1 genome reference.

For each sample, a heat map displaying the proportions' ratio was plotted to visualize the degree of association between mutations located in different amplicons. To simplify the color scale, the highest proportion was placed in the denominator of the ratio, to obtain values in the 0–1 range. Values close to 1 indicate similar proportion of both mutations in the sample.

Statistics

The coefficient of determination (r2) and the Pearson coefficient of correlation (P) were determined to look for any relationship between (i) the average temperature and humidity, (ii) the initial input of the Alpha variant, (iii) the average N1 log10 GC/100,000 inhabitants during the previous five months to the initial detection of the Alpha variant, and the outcompetition rate of the Alpha variant. Additionally, the same coefficients were also assessed between the outcompetition rate of the Alpha variant and the incidence of cases per week in the periods from the first detection of the variant until its dominance (≥95%), and from the first detection till the end of the analysis to include the effect of the deconfinement measures. All statistical analysis were run in SigmaPlot v11.0 (Systat Software Inc, CA, USA).

Results

Emergence and spatiotemporal evolution of the alpha VOC

From November 2020 till April 2021, 14 WWTPs located across the Autonomous Community of Catalonia in Northeast Spain (Table 1 and Fig. 1 A) covering a population of 2,455,668 inhabitants (32% of the total Catalan population) including part of Barcelona and its Metropolitan area, were analyzed for the presence of the Alpha variant.

Fig. 1

Spatiotemporal evolution of the Alpha VOC in the Autonomous Community of Catalonia, Spain. (A) Areas covered by the WWTPs under analysis. (B) First detection was on January 4, 2021, in the PDL WWTP (red), which covers Barcelona and the South Metropolitan area. (C) On January 11, 2021, the variant spread to the nearby areas and was detected in the MRT, RUB, GRA and VDP WWTPs (orange), and to a second focus in Girona where it was detected in GIR and BAY WWTPs (orange). (D) On January 18, 2021, the variant further spread from both focusses into the MIR and FIG WWTPs (yellow), respectively. (E) On January 25, 2021, the spreading continued affecting the SRS, ABR and VEN WWTPs (green). (F) On February 1, 2021, the variant was detected in the BER WWTP (cyan). (G) Finally, on February 8, 2021, the variant was detected in the PUI, WWTP (blue). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

First detection of the Alpha variant was on January 4, 2021, in the WWTP PDL, which covers Barcelona city and the South of its Metropolitan area (Fig. 1B). From this moment on, the variant spread to the nearby areas, and on January 11, 2021, was further detected in the WWTPs MRT (Baix Llobregat), RUB (Vallès Occidental), GRA (Vallès Oriental) and VDP (Alt Penedès), and farther North to the WWTPs GIR (Gironès) and BAY (Pla de l’Estany), representing these latter regions a second focus of virus expansion (Fig. 1C). The next weeks (Fig. 1D and E), the Alpha variant advanced from these two focusses, into the WWTPs MIR (Vallès Occidental) and FIG (Alt Empordà) on January 18, 2021, and further into the WWTPs ABR (Baix Llobregat) and SRS (Vallès Occidental), as well as into the WWTP VEN (Baix Penedès) on January 25, 2021. Finally, on February 1 and 8, 2021, the virus was already detected in the less populous regions covered by the WWTPs BER (Berguedà) and PUI (Cerdanya) in the Pyrenees mountains (Fig. 1F and G). In a period of 6 weeks, the Alpha variant was present in all WWTPs with an average spread of 293 Km2/week, which represents a 17% of the area under study per week, with a maximum of 679 Km2 (39%) in the week of January 11, 2021, and a minimum of 23 Km2 (1%) in the week of February 1, 2021.

S gene markers in the different WWTPs

The presence of the Alpha VOC was confirmed by deep sequencing of the complete S gene in most WWTPs (Table 2 ). Compared to the Wuhan-Hu-1 variant, VOC Alpha is characterized by the occurrence of 9 genetic markers in the S gene (69-70del, 144del, N501Y, A570D, D614G, P681H, T716I, S928A and D118H). Additionally, the replacement E584K has also been associated with some Alpha sequences. We were able to detect a range from 1 to 8 Alpha signature markers in the different WWTPs analyzed (Table 2). The relative abundance of each marker from most to least abundant was D614G > A570D > 69-70del and 144del > D1118H > S982A > P681H and T716I > N501Y > E584K (Table 2).

Table 2

Genetic markers from the S gene detected in the different WWTPs under study.

Temporal dynamics of alpha VOC in the WWTPs

Once the Alpha variant emerged in each WWTP, its temporal evolution was analyzed until its complete dominance was reached. Weekly percentages were estimated using the duplex TaqMan RTqPCR assay for the allelic discrimination described in section 2.4 (Fig. 2 A–K). Independently of the outcompetition rate of the Alpha variant, dominance was achieved in all WWTPs (Fig. 2). On average, 10.73 ± 1.79 weeks were required to reach a ≥95% dominance, with a mode of 12 weeks and a range of 7–14 weeks.

Fig. 2

Temporal evolution of the frequency of the Alpha VOC in each WWTP analyzed. Grey bars represent the relative frequency of the Alpha VOC, measured as the relative presence of the 69–70 deletion, vs. other variants in black, measured as the relative absence of the 69–70 deletion. The red line represents SARS-CoV-2 titers in sewage expressed as Log10 of N1 genome copies per 100,000 inhabitants. The red and green arrows depict the dates of implementation of the confinement and de-confinement measures, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Outcompetition rates for the Alpha variant were estimated through the slope of the regression line between time (in weeks) and its relative abundance (in percentage): the steepest the slope the higher the outcompetition rate (Table 3 and Figure S1). Only those WWTPs with reliable quantitation based on the described criteria in section 2.4 regarding the LOQ were included in the analysis (i.e. PUI, BER and BAY WWTPs were excluded). The mean rate was 6.13 ± 1.39 with a range between 4.12 and 9.78. The slowest rate was observed in PDL, which partially covers Barcelona city and its South Metropolitan area, and the fastest rate in SRS. Remarkably, the first emergence of the Alpha variant was in PDL while SRS was among the last. Although both rates behave as outliers (Figure S2), the outcompetition rate in SRS showed a much higher deviation and was omitted in further analysis.

Table 3

Slope of the Alpha VOC outcompetition curve in each WWTP.

Wastewater treatment plant	Slope of Alpha VOC outcompetition curve
VEN	5.99
ABR	6.28
VDP	6.49
MRT	5.47
RUB	5.98
FIG	6.49
GIR	6.16
GRA	5.49
MIR	5.19
SRS	9.78
PDL	4.12

Despite it has been described that viral loads in upper respiratory samples are 1 Log10 higher for the Alpha variant compared with the B.1.177 variant (Jones et al., 2021), our study did not show a relationship between the increase of the relative abundance of the Alpha variant in sewage and the increase of the total GC per 100,000 inhabitants (Fig. 2). However, the effect of the lockdown, causing a decrease in the overall virus titers, might well be a confounding factor. In contrast, a significant positive correlation (r2 = 0.45; P = 0.67; p = 0.034) was found between the outcompetition rate of the Alpha variant in sewage, and the increase in the weekly number of cases per 100,000 inhabitants, in the period since its initial detection in sewage until its complete dominance (Fig. 3 A). In most WWTPs (VEN, ABR, MRT, FIG, GIR, GRA, MIR and PDL) the deconfinement started before the Alpha variant was dominant. With the aim to include the effect of the deconfinement in all plants, the same correlation analysis was repeated but relating the rate of outcompetition with the incidence increase in the period since the first detection of the Alpha variant until the end of the study (April 12, 2021), when the Alpha variant was dominant in all WWTPs (Fig. 3B). The same trend was observed, with significant positive correlation (r2 = 0.43; P = 0.66; p = 0.038) between both variables.

Fig. 3

Relationship between the Alpha VOC outcompetition rate and the COVID-19 incidence per week. (A) Regression analysis between the outcompetition rate and the incidence per week from the first introduction of the Alpha VOC until its complete dominance. (B) Regression analysis between the outcompetition rate and the incidence per week from the first introduction of the Alpha VOC until the end of the study. The coefficient of determination (r2), the Pearson's coefficient of correlation (P) and its significance (p) are depicted.

The temporal dynamics was also assessed in some WWTPs through the analysis of the genetic markers typical of the Alpha VOC in the S gene (Table 4 ). We selected four WWTPs for their low presence of inhibitors: ABR, VDP, GIR and GRA. In all cases, markers of the Alpha variant were still present at the end of April 2021.

Table 4

Temporal evolution of genetic markers from the S gene in three WWTPs.

Factors influencing the outcompetition rate of alpha VOC

With the aim to explain the different rates of outcompetition of the Alpha variant in sewage (Table 3), we looked for correlations between different variables including environmental factors (i.e. temperature and humidity), the initial input of Alpha VOC, and the mean of the log10 N1 GC/100,000 inhabitants in sewage during the previous five months, as an estimate of population immunity.

No correlations were found either with temperature or humidity (data not shown). In contrast, with the omission of the SRS WWTP, we observed a significant positive correlation (r2 = 0.60; P = 0.77; p = 0.009) between the initial input of the Alpha variant in sewage (i.e., log10 N1 GC/100,000 inhabitants of the Alpha variant on the first detection) and the outcompetition rate (Fig. 4 ).

Fig. 4

Relationship between the SARS-COV-2 levels in sewage at the first day of detection of the Alpha VOC and its outcompetition rate. The coefficient of determination (r2), the Pearson's coefficient of correlation (P) and its significance (p) are depicted.

The mean of the log10 N1 GC/100,000 inhabitants in sewage was figured from August 31, 2020, until the first detection of the Alpha variant in each WWTP (Fig. 5 A) and compared with the total number of clinical cases in this period (Fig. 5B). The relationship between both variables showed a significant positive correlation (r2 = 0.47; P = 0.68; p = 0.02) (Fig. 5C). Next, we used the mean of the log10 N1 GC/100,000 inhabitants for this period, as an indication of previous exposure to SARS-CoV-2 in the population covered by the WWTPs, and we also found a significant positive correlation (r2 = 0.72; P = 0.85; p = 0.002) with the outcompetition rate (Fig. 6 ).

Fig. 5

Relationship between the levels of SARS-CoV-2 genomes in sewage and the numbers of COVID-19 cases in the population. (A) Boxplot distribution of the weekly Log10 N1 GC/L in sewage per 100,000 inhabitants during the period from August 31, 2020, until the first introduction of the Alpha VOC. The 5, 25, 50, 75 and 95 percentiles are represented by the lower whisker, the lower box, the discontinuos line, the upper box and the upper whisker, respectively. Dots correspond to outliers. (B) Total number of cases per 100,000 inhabitants in the same period. (C) Relationship between the total number of cases per 100,0000 inhabitants and the average levels of SARS-CoV-2 N1 GC/100,000 inhabitants in sewage. The coefficient of determination (r2), the Pearson's coefficient of correlation (P) and its significance (p) are depicted.

Fig. 6

Relationship between the SARS-COV-2 levels in sewage during the period from August 31, 2020, until the first introduction of the Alpha VOC, and its outcompetition rate. The coefficient of determination (r2), the Pearson's coefficient of correlation (P) and its significance (p) are depicted.

Dynamics of genetic markers atypical from alpha VOC

Apart from the ten characteristic signatures of Alpha VOC, we found a total of eighteen additional markers distributed in the different WWTPs (Table 2, Table 4). Seven of these mutations (L5F, T73A, S98F, S155I, L216F, D253N and D290Y) were located in the amino terminal domain (NTD) of the spike glycoprotein S, two (S477N and A522S) were located in the receptor binding domain (RBD), and nine (G669V, A706V, K854N, S940A, Q1071L, H1083R, D1084R, G1085R and G1267V) outside relevant regions. Although all these mutations had been previously documented in the GISAID platform (https://www.gisaid.org/) (ElbeBuckland-Merrett, 2017; Shu McCauley, 2017), some relevant features could be anticipated from our data.

Replacements T73A and D253N, which are in the N2 and N5 loops of the NTD, respectively (Cerutti et al., 2021a), were already detected at the beginning of February 2021 and were still present at the end of April in the VDP WWTP (Table 4). These two replacements were likely contained in the same genome molecule as revealed by their association in a proportions' ratio heat map (Figure S3).

Additionally, in the middle of February 2021 at the MIR WWTP, two replacements located in the RBD (S477N and A522S) were identified in one amplicon, meaning that both mutations co-occurred within the same genome molecule. Moreover, these two mutations are likely to be associated with the 69-70del and 144del in the proportion's ratio heat map (Figure S3), suggesting the occurrence of an Alpha variant sequence harboring additional mutations.

Discussion

A clear understanding of the dynamics of spreading of the SARS-CoV-2 is required to design specific mitigation efforts. Few studies have explored how such dynamics may differ in areas with different population densities (Moreno et al., 2020), even though interventions must be addressed taking into consideration time and site characteristics, particularly with the introduction of emerging VOCs.

The first introduction of Alpha VOC in Catalonia, Spain, occurred in late December 2020 (Perramon et al., 2021), and two weeks later confinement measures were implemented. One measure was the prohibition of movements outside municipalities, except traveling to work or to essential appointments, representing a good setting to evaluate the fitness of a VOC and its behavior in a competitive scenario.

It has been extensively proven that genome titers of SARS-CoV-2 in sewage correlate with the number of COVID-19 cases in the population (Chavarria-Miro et al., 2021; Lodderde Roda Husman, 2020; Medema et al., 2020; Peccia et al., 2020). Consequently, information gathered from wastewater surveillance has enabled the study of the spatiotemporal dynamics of Alpha VOC spread in Catalonia from areas with different population densities. Although, we cannot rule out more than one introduction of the Alpha variant in the area under study, confinement measures make this possibility unlikely. Assuming a single introduction of the Alpha variant, through the Barcelona Metropolitan area, it would have spread from the coast towards the interior and to the Pyrenees with an average of one sixth of the available area per week.

Once the Alpha variant was introduced, the rate required to become dominant, i.e., the outcompetition rate, differed between WWTPs. To understand this kinetics, we analyzed the relationship between the outcompetition rate and a series of variables. The differences in the rate of the Alpha variant outcompetition were not explained by environmental variables such as temperature and humidity. At least two factors influencing the virus-host interaction could partially explain the variation in the outcompetition rate: the initial relative abundance of the Alpha variant and the previous circulation of SARS-CoV-2 within the served population. While the potential influence of the initial input is intuitive, i.e., the highest the input of the variant the fastest the rate of outcompetition, the relationship with the previous circulation of SARS-CoV-2 (mainly the B.1.177 variant) within the population is less evident. Although the Alpha variant may escape from antibody responses generated by natural infection, reduction of the neutralization is modest (Supasa et al., 2021). However, the higher the virus circulation the higher the herd immunity, implying a reduction of potential hosts, and a more efficient transmissibility of the Alpha variant (Davies et al., 2021; Volz et al., 2021) might represent a clear advantage.

Additionally, the co-circulation of the Alpha variant with other variants able to efficiently escape neutralization by antibodies generated during the natural infection or by vaccination (the average of two-dose vaccination at the end of the study was of 5% ± 2% of the total population), or with high transmissibility (i.e., Beta, Gamma and Delta) could also influence the spread kinetics of the Alpha variant. During the period under study the Beta and Gamma VOCs were identified in 2% and 4% of the sequenced clinical samples in Catalonia, and were mostly detected in populations served by the PDL, MIR, SRS, RUB, and GRA WWTPs, (https://canalsalut.gencat.cat/web/.content/_A-Z/C/coronavirus-2019-ncov/material-divulgatiu/informe-vigilancia-variants-genomiques-sarscov-2.pdf), which were the ones showing the lower outcompetition rates, with the exception of SRS.

The extremely fast outcompetition rate of the Alpha variant in the SRS WWTP remains unexplained. It could respond to a second massive introduction of the variant between six and seven weeks after the original one (Fig. 2J). The outstanding fast outcompetition rate in this case, did not correlate with a great increase in the number of cases (data not shown) but it could be supported by a surge of asymptomatic infections. Furthermore, the SRS WWTP covers a densely populated region with high industrial and commercial activity, which opens the possibility of employees shedding the virus at their working place (e.g. during the incubation period) but with case reporting at their hometowns covered by other WWTPs. Additionally, differences in the composite sample collection, microbial community, chemical pollution and pH of sewage, residence time of sewage in the pipe network, etc., between the WWTPs, may also influence the level of GC/L and contribute to certain degree of bias.

Deep-sequencing analysis of sewage may provide additional information on the emergence and evolution of VOCs and/or mutations of concern (Fontenele et al., 2021; Izquierdo-Lara et al., 2021). Mutations T73A and D253N are located in loops N2 and N5, respectively, of the NTD of the S glycoprotein of SARS-CoV-2 (Cerutti et al., 2021a). While loop N5 is integrated in the antigenic supersite, loop N2 is outside. However, it has been suggested that mutations in the N2 loop could significantly impact the conformation of the NTD supersite (Cerutti et al., 2021a). In consequence, the occurrence of both mutations in the same genome in two samples from the VDP WWTP collected 11 weeks apart, may suggest a certain advantage of the virus harboring them. Similarly, mutations S477N and A522S which were present in the same genome molecule, could represent a combination of concern. In the period under study, replacement S477N was one of the most common in the GISAID database, while A522S although frequent was much less common. In fact, both mutations had been described to be among the top 12 mutation positions in the RBD (Cerutti et al., 2021b; Wang et al., 2021). Mutation S477N has been described to induce resistance to neutralization by multiple monoclonal antibodies and, in the presence of the S514F replacement, to human immune sera as well (Liu et al., 2021). The effect of the combination of the S477N and A522S mutations on the resistance phenotype to immune sera has not yet been elucidated, but it cannot be rule out. However, these mutations were only detected in mid-February, suggesting that they would not provide any advantage, or even they could be associated with a loss of fitness, in the absence of herd immunity. The identification of VOCs in sewage is not as easy as it can be anticipated due to many difficulties such as: i) occurrence of many mutation signatures shared between variants, ii) lack of robust algorithms for the identification of the occurrence of distant mutations and/or deletions in a genome and iii) technical problems such as the lack of amplicon coverage homogenization, particularly when virus titers are low. Despite these technical difficulties, the identification of signatures (mutations and indels) of concern and the study of their spatiotemporal evolution in sewage could be highly valuable for public health management.

Nevertheless, in spite of its intrinsic technical limitations, Wastewater Based Epidemiology offers a valuable overview of the dynamics of virus circulation among the population in a given area.

Conclusions

Despite the usefulness of clinical-based analysis to mitigate SARS-CoV-2 spread, asymptomatic infections are mostly undetected. In consequence, wastewater-based epidemiology becomes a good complementary tool for a broad virus surveillance including variants characterization. This study provides the evolution of viral load titers in 14 wastewater treatment plants serving populations with different demographic features, as well as the spaciotemporal dynamics of the Alpha VOC spread. In a period of 6 weeks the Alpha VOC was present in all the studied WWTPs, and became dominant in 11 weeks on average. However, the pace of spread was not homogeneous across the WWTPs. Obviously, the pace was faster when the Alpha VOC introduction in a population was stronger, but also when the previous circulation of SARS-CoV-2 was higher. This latter point could be related with the advantage of a more transmissible variant in a population with a certain level of herd immunity. The aproximation used for the estimation of the Alpha VOC outcompetition rate could be applied to other VOCs, broadening the WBE-derived results which in consequence may further contribute to the understanding of the SARS-CoV-2 epidemics.

When the coverage of clinical samples consensus sequencing is insufficient, deep-sequencing of sewage samples, albeit its limitations, is an excellent alternative. In this study, we have identified the co-occurrence of two pairs of mutations (T73 A/D253N, and S477N/A522S) which could be of concern, reinforcing the importance of the WBE for the early detection of potential novel variants.

Funding sources

This work was partially funded by the Catalan Agency for Water (), the Catalan (ASPCAT) from the , and the Health Innovation Program from the General Research Directorate (DGRIS) of the , and the Spanish Centre for Technological and Industrial Development () from the Spanish and Business (grant IDI-20200297).

Competing financial interests

The authors declare they have no actual or potential competing financial interests.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.