Wastewater sampling

From January to May 2023, weekly wastewater samples were collected from aircraft, Arlanda Airport Terminal 5, Måby station (covering the airport region), Käppala WWTP, Bromma WWTP and Henriksdal WWTP (Figs. S1 and S2, supplementary material). Stockholm has an urban population of around 1.6 million and a metropolitan population of ~2.4 million. Samples were collected from three main municipal wastewater treatment plants in Stockholm: Bromma WWTP (serving around 377,500 residents, 18% of the population), Henriksdal WWTP (serving about 862,100 residents, 41% of the population) and Käppala WWTP (serving around 700,000 residents, 33% of the population). Bromma WWTP has three inlets: Hässelby, Riksby and Järva, while Henriksdal WWTP has two inlets: Sickla and Henriksdal. All airport wastewater goes to Käppala (Fig. S1) and this WWTP has only one inlet. The samples from the six inlets and the airport were flow-compensated and ~500 mL of each sample was collected over 24 h from Monday to Tuesday every week, using stationary flow-proportional samplers SP5 B (MAXX Mess- und Probenahmetechnik GmbH, Germany) in Henriksdal WWTP, TP5 W (MAXX Mess- und Probenahmetechnik GmbH, Germany) in Bromma WWTP and Efconomy (Efcon Water B.V., The Netherlands) in Käppala WWTP.

Stockholm’s samples were collected from week 1 to 22 of 2023 while aircraft samples from week 2 to 18, except for the 3rd and 11th weeks when flights were cancelled. A direct flight from China carrying 312 passengers and with a duration of ~9 h was weekly monitored, on Fridays from week 2 to 9 and on Tuesdays from week 10 to 18. Aircraft wastewater samples were collected after thorough mixing during transport from the aircraft to the lavatory service truck. Samples were then taken from the top of the lavatory service truck immediately after the transfer. Grab samples, reflecting the microbiological status at the time of collection, were taken from Stockholm’s Arlanda Terminal 5 with data available in supplementary material (Figs. S4 and S5). All samples were transported to the lab on cooling boxes and with ice packages. Samples from aircraft and airport (Måby and Arlanda terminal 5) were processed for concentration and RNA extraction in the same day (less than 24 h after sampling), while samples from the six inlets of the WWTP were processed within 24 h after sampling.

Wastewater concentration and RNA extraction

Wastewater samples were concentrated using a Maxwell RSC Enviro TNA Promega Kit following the manufacturer’s instructions, with some exceptions. From the 500 ml of wastewater sample collected weekly, 40 ml aliquots were used for analysis and treated with a protease solution and centrifugated to remove precipitated proteins and solids present in the sample. Then, the supernatant was filtrated and eluted to 500 µL using a column-based system and then loaded into a cartridge provided by the kit. Afterwards, total nucleic acid (TNA) was extracted using Maxwell RSC Instrument (Promega Biotech AB, Sweden) and Maxwell RSC Pure Food GMO programme was selected. The elution volume was 80 µL using nuclease-free water. Two independent biological replicates were concentrated and RNA extracted per sample. Two tap water samples were used as negative controls in each extraction set (extraction of 16 samples). To ensure unbiased analysis, samples were labelled with specific codes for each location and week, and these codes were randomised during the wastewater concentration and TNA extraction processes.

Due to the relatively low SARS-CoV-2 content in aircraft samples, a spatial composite of aircraft wastewater was implemented from week 5 onwards to enhance the likelihood of positive detection and increasing the amount of virus before RT-qPCR analysis. Four independent wastewater samples, each with a volume of 40 mL, were concentrated and eluted to a final volume of 500 µL using the Promega kit. The four resulting concentrates of 500 µL each were then divided into two groups, and the samples within each group were combined, yielding ~1000 µL per group. These combined samples were subsequently used for RNA extraction, producing two independent biological replicates, each with a volume of 80 µL.

Weekly monitoring of SARS-CoV-2

SARS-CoV-2 contents were determined via reverse transcriptase quantitative polymerase chain reaction (RT-qPCR)^35,36. The reaction was performed using SYBR Green one-step kit (Bio-Rad) according to the manufacturer’s instructions, with the modification of adding 2 μL of 4 mg/ml Bovine Serum Albumin (BSA) (Thermo Scientific™) to reduce PCR inhibitors and enhance efficacy, for a final reaction volume of 20 μL. For SARS-CoV-2 detection, N3 primers targeting the nucleocapsid (N) protein were used, FW: 5′-GGGAGCCTTGAATACACCAAAA-3′ and RV:5′-TGTAGCACGATTGCAGCATTG-3′³⁷. Pepper mild mottle virus (PMMoV) was also quantified in the samples and used for data normalisation as previously described³⁶. The primers used for PMMoV quantification were forward, 5′-GAGTGGTTTGACCTTAACGTTTGA-3′; reverse, 5′-TTGTCGGTTGCAATGCAAGT-3′³⁸. Nuclease-free water was included as no-template control (NTC) for all qPCR reactions. Two tap water samples were concentrated, and RNA was extracted and analysed alongside the weekly samples to assess potential contamination during handling. SARS-CoV-2 DNA (2019-nCoV_N_Positive Control, IDT, Cat. 10006625), and a constructed plasmid containing the appropriate target for PMMoV (IDT, custom MiniGene 25-500 bp) were used as positive controls and to create the standard curves. A control sample of RNA from wastewater (cross-plate controls) with known SARS-CoV-2 and PMMoV concentrations was used in each qPCR analysis for reproducibility and quality control. Cross-plate controls were stored at −80 °C in aliquots. A new batch was prepared when the quantification cycle (Cq) shift of 0.5–1 was detected, with each batch typically used for up to 8 weeks. The cross-plate control was employed to evaluate the long-term precision, which refers to the variation in results between runs, in accordance with MIQE guidelines³⁹. The standard deviation of the cross-plate controls was calculated between runs. Results were accepted if the standard deviation of the Cq values was less than 0.5. Two independent technical replicates were measured for each biological replicate and each control sample (positive, negative and cross-plate controls) using qPCR. Results were accepted if the standard deviation of Cq values between technical replicates was less than 0.5. Thermal cycling (50 °C for 10 min, 95 °C for 30 s, followed by 45 cycles of 95 °C for 10 s and 60 °C for 30 s) and melting curve detection were performed (65 °C to 95 °C with an increment of 0.5 °C for 5 s) on CFX96 Touch System (Bio-Rad). Reactions were considered positive if the fluorescence crossed the established threshold before 40 cycles (if Ct was less than 40) and if a single melting peak was observed at the correct temperature. The threshold was set automatically using CFX Manager™ Software. The standard curves yielded calculated efficiencies of 98% and 90% for the N3 and PMMoV qPCR assays, respectively, with corresponding slopes of −3.38 and −3.59. The coefficient of determination (R²) values were 0.999 and 0.998, respectively. The LOD and limit of quantification (LOQ) for the qPCR assays was determined using the standard deviation method, which involves the standard deviation of the response (y-intercepts of the regression lines) and the slope of the calibration curve. The LOD and LOQ of the qPCR assays used for N3 were 0.4 copies/ml of wastewater sample and 5.3 copies/ml, respectively, and for PMMoV, they were 0.4 copies/ml and 4.1 copies/ml, respectively. Inhibition testing was conducted utilising the Cq dilution method, as recommended by the MIQE guidelines³⁹. The results demonstrated high efficiencies and strong R² correlation coefficients for the standard curve. An example of the inhibition test is provided in the supplementary material (Fig. S9). A sample was considered positive if at least three out of four measurements were positive, including two technical replicates for each of the two biological replicates. If only two out of four measurements were positive, the samples were re-tested by qPCR. If the re-test yielded the same result and the melting curves were inconclusive, the samples were concentrated and tested again by qPCR.

SARS-CoV-2 concentrations in aircraft samples were near the detection limit for several weeks during the monitoring period. For this reason, we evaluated whether the detection limit could be improved by using TaqMan method or if there was any difference among the well-known primers targeting the N protein (N1, N2 and N3 primes). The methodology and results are presented in supplementary material.

Calculations

The SARS-CoV-2 contents in this study are presented as either N-gene copy number adjusted for variations in PMMoV contents per week (PMMoV factor, Fig. 2) or N-gene copies per PMMoV gene copies ×10⁴ (N-gene/PMMoV ratio, Fig. 1). The PMMoV factor and N-gene/PMMoV ratio calculation methods were previously described by Perez-Zabaleta et al.³⁶. The PMMoV factor calculation method, which adjusts for changes in flow rates, has been utilised to analyse data from Käppala and Stockholm. This approach combines data from various inlets to provide the total SARS-CoV-2 detected in the Stockholm region. The N-gene/PMMoV method was not used to plot Stockholm data due to its inability to account for variations in flow rates at each testing location, which could lead to inaccurate results. For locations such as aircraft, Arlanda Terminal 5 and Måby (airport), where flow rate data is unavailable, the N-gene/PMMoV ratio method was used to present the results. However, these results were kept separate and not combined with each other. For reference, Fig. S6 was included in the supplementary material to present the data of the five locations without normalisation (PMMoV or flow rate).

DNA sequencing, variant calling and lineage determination

After the wastewater concentration and RNA extraction steps, 25 µL of purified RNA was shipped to the Uppsala Genome Center (Science for Life Laboratory, Dept. of Immunology, Genetics and Pathology, Uppsala University) for sequencing. The two biological replicates were mixed in equal amounts (12.5 µL) and sent for sequencing per sampling point. These sampling points included the six inlets of the WWTPs (Henriksdal, Sickla, Hässelby, Järva, Riksby and Käppala), aircraft, Arlanda Terminal 5 and Måby (airport region) samples. For Stockholm composite samples, 4.5 µL of each WWTP inlet were pooled, resulting in a volume of 27 µL of sample per week. For weeks 7, 8 and 9, the mixed sample from Stockholm was not sequenced, and data was pooled in silico. Additionally, in silico pools including all six inlets or five inlets excluding Käppala, were created when needed and specified in the results.

The samples were then reverse-transcribed with SuperScriptVILO cDNA synthesis (Thermo Fisher Scientific, Waltham, MA, USA) and the libraries were prepared using the AmpliSeq SARS-CoV-2 panel (Thermo Fisher Scientific) on two S5 540 chips. The libraries were barcoded and pooled into 32-plex and sequenced on the Ion S5XL system (Thermo Fisher Scientific). Processing raw sequencing reads, and single nucleotide variant (SNV) calling was done through the Torrent Suit AmpliSeq SARS-CoV-2 pipeline according to the manufacturer’s instructions. The resulting BAM files were used for post-processing with the Freyja pipeline, v.1.4.8⁸. Briefly, ‘Freyja variants’ were used to track SNVs against the Wuhan-2019-nCoV reference genome. Then, to convert the depth of each SNV into a likely relative abundance of PANGO lineages⁴⁰, ‘Freyja demix’ was used with barcode version 12_12_2023-00-48, a depth cutoff of 10x, a minimum relative abundance (eps) of 0.001 and only confirmed variants^41,42,43. Sequences are deposited at ENA under https://www.ebi.ac.uk/ena/browser/view/PRJEB61810. For statistical analyses, only variants that could be fully assigned by Freyja, without the ambiguity markers ‘Misc’ or ‘-like’ were considered. Figures 3, 4 and S7, S8 were created in R v. 4.3.1, with libraries VennDiagram 1.7.3 and UpSetR 1.4.0.

Clinical SARS-CoV-2 variant data

Clinical infection data on COVID-19 patients were only available from Stockholm city since no epidemiological data corresponding to the aircraft, airport and Käppala were available. Test recommendations changed in February 2023, from including patients and healthcare personnel to patients only. Sequences from clinical cases are available in the Source Data file and GISAID EPI_SET_250516ap, where a proportion of positive samples from patients were sequenced as part of the genomic surveillance programme.

Statistical analyses

Pearson correlations were performed to determine if there was any correspondence in SARS-CoV-2 content between Måby and Käppala (Fig. 1, N-gene/PMMoV ratio) or Käppala and Stockholm (Fig. 2, normalised SARS-CoV-2 weekly viral load). P < 0.05 were considered statistically significant. All the statistical analyses were conducted in Prism version 10 (GraphPad Software, CA, USA). Standard deviations (s) were calculated by Eq. (1) and the variance (s²) using Eq. (2), where: (x_i) represents each data point, (µ) is the mean of the data and (N) is the number of data points.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.