Sample collection

This study includes data from samples collected over the period of January 1st 2022 to March 31st 2023. WW samples from two airport terminals at Toronto Pearson (Terminal 1 and Terminal 3) were collected using passive torpedo samplers (soaked for 24 h each), until March 21st 2022 after which Terminal samples were collected using autosamplers that took one grab sample each half hour over a 24 h period to create a composite sample. Pooled aircraft sewage was also first collected from January 1st 2022 until October 21st 2022, using passive torpedo samplers. After October 21st 2022, pooled aircraft sewage was collected using an autosampler in the same way as the Terminal samples yielding 24 h composite samples. The switch to autosamplers was reliant on the installation of autosamplers at each location and the terminals were prioritized due to the WW flow at these sites damaging the passive samplers. Samples from municipal WWTPs were collected similarly by autosampler in the case of sites 1–6 and 9 or by grab sampling (single timepoint collection) in the case of sites 7 and 8, at pumping stations or pipes. Sample metadata including sample type, coverage statistics and sample preparation method is available in Supplementary Table ^S1.

Library preparation and WGS

As part of the Ontario WSI, for each site included in this study, WGS was conducted by one of three partners at University of Guelph, University of Waterloo, and University of Western Ontario. The following section details which samples were processed by each lab, and how the procedures and reagents used differed.

WW samples were subjected to RT-qPCR to measure the quantity of SARS-CoV-2 present in the samples prior to sending them to one of the sequencing labs, or for the Terminal and Aircraft Sewage samples, the RT-qPCR was performed directly in the sequencing lab. All sequenced samples are summarized in Table 1. The SARS-CoV-2 concentrations in airport samples were stable throughout the sampling period (Ct scores between 29 and 31, Supplementary Table ^S1). During this study, municipal samples that were part of the WSI program with Ct scores above 35 were not sequenced, thus all municipal samples included here have Ct scores≤35. Gene copy measurements for the entire Ontario WW qPCR dataset were published as copies/L values for all municipal samples included in this study³⁸.

Table 1

Summary of samples.

Site Name	Region	# samples passing QC (Ct<35, BOC10>20%)	# samples collected/week	Sample type
Pooled Aircraft	Airport	130	5	Passive Torpedo Sampler 24 h
		91		Autosampler 24 h
Terminal 1		44	5	Passive Torpedo Sampler 24 h
		233		Autosampler 24 h
Terminal 3		46	5	Passive Torpedo Sampler 24 h
		246		Autosampler 24 h
Leslie Street PS	York	55	1	Grab
Warden 407		53	1	Grab
Peel OCF/Humber AMF		62	1	Autosampler 24 h
Clarkson	Peel	60	1	Autosampler 24 h
G.E. Booth		61	1	Autosampler 24 h
Humber (Etobicoke)	Toronto	76	1	Autosampler 24 h
Ashbridges Bay		56	1	Autosampler 24 h
North Toronto TP		54	1	Autosampler 24 h
Highland Creek		79	1	Autosampler 24 h

Open in a separate window

Table summarizing all samples analyzed for this study. Site names and regions have been shortened. Region names will be used to define groups of samples as shown throughout results. Number of samples sequenced minus samples that were discarded because the breadth of coverage (BOC) at 10 reads deep or greater was less than 20% is counted in column “# samples passing QC (Ct<35, BOC10>20%)”. Sample type refers to how the sewage was collected.

One sample per weekday (Monday - Friday excluding statutory holidays) was collected from Toronto Pearson Terminal 1 and Terminal 3 for sequencing at University of Waterloo. One sample per week was collected for sequencing at University of Waterloo from each of the following municipal WWTPs, pumping stations, or pipes: Leslie Street Pumping Station - York Sewershed (Site 7), Warden and 407 - York Sewershed (Site 8), Humber AMF Pumping Station - York Sewershed (Site 9), Clarkson South-Peel Water Pollution Control Plant (Site 1), and G. E. Booth (Lakeview) WWTP (Site 2).

In preparation for sequencing, the viral content of 20 mL of a WW sample was concentrated using Nanotrap Microbiome A Particles (Ceres Nanosciences) following the manufacturer’s protocol “Manual Nanotrap^® Wastewater Protocol using QIAGEN QIAamp^® Viral RNA Mini Kit”, section A. “Viral Capture” with the following modifications: whenever a magnetic stand was called for, the sample with beads was instead centrifuged at 8,000 × g at 4 °C for 4 min, and the pellet was resuspended in a small volume of remaining supernatant instead of water to transfer the sample to a smaller tube. In the case of Terminal 1 and 3 samples prior to March 21st 2022 one of 3 gauzes from the torpedo passive sampler was first shaken in 40 mL sterile PBS with 0.05% Tween 80 and antifoam reagent, wrung out, and the PBS with resuspended biosolids was used as the input material for the Nanotrap Microbiome A Particles (Ceres Nanosciences) concentration. The beads were resuspended in 700 µL of prepared Lysis Buffer (RLT+2ME) for RNA extraction automated on a QIAcube Connect using a RNeasy mini kit (Qiagen ID: 74116). RNA was then reverse transcribed using LunaScript RT SuperMix Kit (NEB #E3010) and amplified using Q5^® High-Fidelity 2× Master Mix (NEB cat# M0492L) and the ARTIC V4.1 NCOV-2019 Panel of primers (IDT; 10008554, designed by the ARTIC network) resulting in ~400 bp amplicons. Amplicon libraries were prepared using the Illumina DNA Library Prep kit (20060059) with an additional 0.8× AMPure XP (Beckman Coulter) single sided bead cleanup prior to tagmentation. Paired-end (2×250 bp) sequencing of the libraries was performed using a MiSeq (Illumina).

One sample per day Monday-Friday (excluding statutory holidays) was collected from the triturator for sequencing at University of Guelph. Pooled aircraft sewage samples were prepared and sequenced as described previously³⁹. Briefly, samples from passive samplers were preprocessed by submerging one of the embedded gauzes in 50 mL sterile PBS with 0.05% Tween 80 and antifoam reagent and stomached for 2 min. Viral particles in the filtrate from passive sampler and in the liquid samples from autosamplers were concentrated using Nanotrap Microbiome A particles (Ceres Nanosciences). RNA extraction was performed using the QIAamp Viral RNA extraction Mini kit (QIAGEN) carried out on the QIAcube (QIAGEN) instrument according to the manufacturer’s extraction method. Complementary DNA was generated using the SuperScript™ IV First-Strand Synthesis System (Thermo Fisher) and amplified using the same ARTIC V4.1 primers and protocols as above for the University of Waterloo. The amplicon libraries were prepared using the Nextera XT DNA library prep kit (Illumina). Paired end (2×150 bp) sequencing of the libraries was performed on the Illumina MiniSeq system. A maximum of 24 samples were loaded per run.

One sample per week was sequenced at the University of Western Ontario from each of the four Toronto WWTPs: Humber (Etobicoke) Water Pollution Control Plant (Site 3), Main (Toronto-Ashbridges Bay) Water Pollution Control Plant (Site 4), North Toronto Treatment Plant (Site 5), and Highland Creek (Scarborough) Water Pollution Control Plant (Site 6). Ceres beads were not used to concentrate the virus, instead a 30 mL aliquot of the composite influent WW sample was centrifuged at 4200 × g at 4 °C for 20 min. The supernatant was discarded, except for approximately 500 µl used to re-suspend the solids pellet. The re-suspended solids pellet was used for RNA extraction as described above for University of Waterloo. CDNA, PCR amplification and Library preparation was performed as described above for University of Guelph. Paired-end (2×300 bp) sequencing was performed using an Illumina MiSeq.

Data processing and frequency prediction

Fastq files for all samples were pre-processed using cutadapt to remove adapter sequences, minimap2 to map reads to a SARS-CoV-2 reference genome (NC_045512.2 – Wuhan-Hu-1), and samtools to generate coverage and mutation frequency data at all sites along the genome^40–42. Mapped genomes with breadth of coverage (BOC)>20% at a minimum read depth of 10 were then processed using Alcov, a tool which uses unique and shared mutations to predict frequencies of SARS-CoV-2 lineages³¹. Shared mutations were used to predict presence of sets of lineages while unique mutations were used to differentiate between individual sub-lineages. Mutations were only considered if the read depth at that location was greater than 10. Constellations defining each lineage in Alcov were built to include all mutations differing from the reference genome which occur in ≥90% of the specific lineage sequences on covSPECTRUM³³. Alcov was used to generate heatmaps and csv tables of predicted frequencies of lineages as well as mutation frequencies, from which the values presented in this article were extracted and plotted.

A benchmarking study has been published comparing Alcov to eight other lineage prediction tools designed for use with WW sequencing⁴³. This study concluded that Alcov performed well compared to these tools including well-known tools such as Freyja⁴³.

Airport WW can provide an early warning of emerging VOCs

The raw sequencing data has been deposited to the Sequence Read Archive (SRA) under BioProject PRJNA1088471. Supplementary Table ^S2 shows processed lineage frequency results that were used to generate plots.

Emerging lineages and mutations of interest were tracked to compare first instances in airport WW, surrounding municipal WW, and clinical data where available. The lineages included in this study were BQ*, BA.2.75*, BF*, XBB*, XBB.1.5*, XBB.1.16* and XBB.1.9*. Only the sequencing results passing QC metrics as defined in the methods section were included in the following plots. A frequency of zero in these plots does not mean there was no SARS-CoV-2 detected in these samples, but rather that the specific lineage being queried was not found in the sample. Summed detection data for each of these lineages was plotted in Supplemental Figures ^S1 through S7 against epidemiological week (epi week), and the first detections were compiled and compared (Tables 2 and ​and3).3). Summed detection data was plotted in Fig. 2 against an X-axis normalizing all timelines of lineage emergence in terms of the first clinical detection as reported by PHO so that the WW detections could be visualized by lead time^27,46–51.

When all BQ* lineage predictions were compiled for all the sites (Fig. 2C, Figure ^S1) results showed that the first detections occurred in the airport terminals in epi weeks 22 and 33, 15 and 4 weeks prior to PHO reported clinical detection, though the first detection was at a frequency less than 5%. Then at epi week 34 the pooled aircraft sewage had a small signal which continued to be the main detection until the Ontario clinical data first reported BQ* in epi week 37⁴⁷. The surrounding municipal sites did not have a detection of BQ* until a small detection in epi week 40 (Toronto) and a slightly higher detection in epi week 41 (York region), both after the first PHO clinical sequence was collected.

BA.2.75* emerged in the Ontario clinical data in epi week 29 (Fig. 2A, Figure ^S2)⁴⁶. Airport terminal BA.2.75* detections in WW occurred in epi weeks 17 and 28, 12 weeks and 1 week prior to PHO reported clinical detection, and the first detection in surrounding municipal samples was in York in epi week 23. All of these detections prior to the clinical first detections were at frequencies less than 5%. During post-processing of the data with Alcov there were some low frequency detections of BA.2.75 as early as February 14th 2022 in York region municipal WW and Terminal WW samples but these were not likely to be true detections. The coverage of samples with low frequency BA.2.75 detections in February and March 2022 were all below 60%, and these samples were predicted to contain parent lineage BA.2 when they were initially reported to PHUs and Toronto Pearson. BA.2.75 was not designated as a Pango lineage until June 2022, with detections earlier than June 2nd 2022 only reported from India⁵². WGS from India early in the pandemic was infrequently submitted and often had ambiguous sequence calls, causing scientists to recommend excluding early sequences from India when analyzing GISAID WGS data^53,54. It is more likely that the low frequency detections that occurred in our dataset were another closely related BA.2* lineage so we have analyzed only detections which occur within 16 weeks of the first clinical detection by PHO through the end of the study period.

BF* lineages were detected in the airport terminals in epi weeks 7 and 10, 16 and 13 weeks ahead of PHO clinical detection respectively, and then in Peel (epi week 12), York (epi-week 21), in Toronto in epi week 21 and multiple times in the airport terminals between epi weeks 21 and 24 before Ontario’s first clinical case was reported in epi week 23 (Fig. 2B, Figure ^S3)²⁷.

Open in a separate window

Fig. 2

Frequency of lineages of interest from the de-mixing performed by Alcov (log sale). Horizontal dashed lines represent the 5% and 1% detection thresholds. Each point represents a single WW sample. Points with a frequency less than 1% are dimmed, as these are too low to consider a lineage as “detected”. The x-axis represents the weeks since the first PHO reported Ontario clinical case, with negative values representing detections in WW prior to clinical detections. Each panel represents a single lineage from the following list: BA.2.75* (A), BF* (B), BQ* (C), XBB.1.16* (D), XBB.1.5* (E), XBB.1.9* (F), XBB* (G).

When XBB* lineages were analyzed as a group, detection occurred in Terminal and pooled aircraft WW during epi week 39, two weeks prior to the PHO initial clinical detection, and then in municipal WW first in Toronto in epi week 45, though consistent detection did not begin in municipal WW until epi week 47. The first clinical detection in Ontario took place during epi week 41 (Fig. 2G, Figure S4)⁴⁷.

In addition to analysis of the XBB* lineages as a group, three specific sub-lineages of XBB were investigated: XBB.1.5*, XBB.1.16*, and XBB.1.9* as these seemed to be the most persistent. The first detection of an XBB.1.5* sub-lineage was in an airport terminal sample in epi week 39, 10 weeks prior to PHO initial clinical detection, and 8 weeks before the first municipal detection in a Toronto sample in epi week 47. There was another detection of XBB.1.5* in epi week 49, and the first clinical detection, also in epi week 49 ⁴⁹. Consistent detection of XBB.1.5* in municipal and airport WW samples was not observed until epi week 1 of 2023 (Fig. 2E, Figure S5).

The first detections of XBB.1.16* were in pooled aircraft sewage during epi week 9, followed by detection in the airport terminals and the first clinical detection in epi week 10⁵⁰, . The first detection in municipal WW occurred in epi week 11 in Toronto (Fig. 2D, Figure S6).

For XBB.1.9*, the first detections were all in the pooled aircraft sewage samples from epi weeks 2–6. The first clinical detection in Ontario did not take place until epi week 5, 3 weeks after the first airport WW detection, whereas the first detection in municipal WW was not until epi week 9 in Peel Region (Fig. 2F, Figure S7)⁵¹.

Cumulative detections above 1% and 5% over time per WW site were plotted in Fig. 3 to visualize the signal consistency after detection of each lineage (Fig. 3A-G). A stepwise increase in the cumulative detections each time a sample was taken is expected if a lineage was introduced into a community and remained present, either by spreading to other individuals or continuing to be shed by the original individuals, as each sample would have a detection of that lineage. More than one sample was collected each week for all sample groups represented in Fig. 3 so consistent detection appears as a line which increases on the y-axis (detections) faster than on the x-axis (weeks). Examples of consistent detection after initial detection were seen in the case of Toronto (purple) BF* (Fig. 3B) and York (pink) BF*, BQ* and XBB.1.5* (Fig. 3B, C and E). In most cases airport terminals 1 and 3 (teal) and pooled aircraft sewage (orange) detections plateaued after initial detection which means the initial detections were isolated instances and were not followed by consistent detection of the same lineage in the following weeks. After the clinical detection the airport WW sites showed consistent detection of the lineages.

Open in a separate window

Fig. 3

Cumulative detection probabilities for each lineage and location. Dotted lines represent detections with a 0.01 (1%) detection threshold, and solid lines represent detections with a threshold of 0.05 (5%). The x-axis represents the weeks since the first PHO reported Ontario clinical case, with negative values representing detections in WW prior to clinical detections. Each panel represents a single lineage from the following list: BA.2.75* (A), BF* (B), BQ* (C), XBB.1.16* (D), XBB.1.5* (E), XBB.1.9* (F), XBB* (G).

GISAID EpiCoV™ is a searchable database for SARS-CoV-2 sequences that is accessible only to registered users. It was possible in June 2024 to filter the Ontario clinical sequences in this database by a sub-set of the lineages of interest focused on in this study, and the found first detections are summarized in Table 3²⁸. PHO-submitted clinical data on GISAID included detections prior to those published in the risk assessment reports. When compared to the GISAID data that was available, Airport WW did not provide as much lead time.

Clinical first detections from PHO risk assessment reports were used as the primary dataset for comparison but clinical first detections from GISAID were also considered if available²⁸. Sequence submissions which were available at the time of publication as representatives of a lineage on GISAID may not have been classified as that lineage during the study period when these lineages were emerging.

The country in which a lineage of interest originated could inform the possible routes of transmission before arrival in Ontario. Countries of origin were collected from the pango-designation GitHub and are as follows: BQ.1 originated in Nigeria, BA.2.75 was initially mainly found in India, BF.1 was initially found in England, Denmark, Spain and Scotland, XBB was initially found in USA and Singapore, XBB.1.5 originated in the USA, XBB.1.16 was initially found in India, USA, Singapore, and Europe, and XBB.1.9.1 was initially found in Indonesia, Singapore, Malaysia and England^44,45. XBB*, XBB.1.5* and XBB.1.16* were associated with or originated in the USA, while other lineages were associated with or originated from countries on continents other than North America.

In Table 4 first detections were summarized statistically. Values were summarized by location, across all lineages. Frequencies greater than 1% or 5% were each used as thresholds to consider a lineage detected in the analysis. The airport terminals had the earliest average first detection lead time, with airplane sewage having the second earliest average lead time, but with lower variance and thus aircraft sewage had more consistent detections.

The lead time of each site compared to the PHO first detect or GISAID first detect was summarized statistically in Table 4, considering only the lineages which were searchable on GISAID for that portion of the table (Table 3)²⁸. The values in this table are represented in days for the GISAID comparisons and weeks for the PHO report comparisons. While there were still some instances where WW detections preceded the GISAID detections which can be seen in the earliest first detection column, the mean first detections were for the most part positive numbers meaning that the WW detections usually occurred after the first sample on GISAID for a lineage.

Clinical sequences available on GISAID represent earlier cases of emerging variants

Throughout the study the clinical first detections as found in PHO reports were primarily used to measure lead time. The GISAID clinical sequences submitted by PHO may not reflect data availability at the time that these lineages were emerging as the data may have been transformed or re-analyzed since that time, while the PHO reports were accessed as they are referenced. It is also unlikely that end users of this data would search for emerging lineages on their own using the GISAID database. A study on variant detection by WGS in airport WW performed in Singapore reported even longer lead times for Omicron sub-variants and similarly used only local public health records for comparison²³. For these reasons, the usefulness of WW surveillance at transportation hubs for informing public health of emerging lineages is not decreased by the shorter (or lack of) lead times when compared to GISAID data.

The presence of earlier GISAID clinical datapoints for the lineages BA.2.75*, XBB*, XBB.1.9* and XBB.1.16* indicated that the lineages were present in low abundance in Ontario prior to detection at Toronto Pearson (Tables 3 and ​and4).4). In cases where the countries of origin included the USA, it is likely that cross-border transmission occurred at land border crossings as often as or more frequently than through the international airport, which would provide a possible explanation. BA.2.75* was an interesting case however with only one very low frequency detection occurring before the GISAID detection, because it originated in India, so had potential to arrive in Ontario from overseas carried by airline travelers.

Transient nature of airport WW and accuracy of first detections

Some considerations should be made for the transient populations contributing to the airport sampling sites. Terminal samples provided information on both the transient population of travelers (people who may have been terminating their journey in Southern Ontario or connecting onward to other destinations) and the more consistent population of airport staff and associates who worked in close proximity to the travelers. Pooled aircraft sewage on the other hand mostly represented the transient populations made up of travelers and mobile aircraft staff. For both airport sampling sites, multiple signal spikes occurred which dropped again before a new lineage gained a foothold in Ontario and more consistent detection of a lineage in municipal WW occurred. Each spike could represent a single case or a few cases passing through the airport, and there is no way to know the ultimate destination of the person or persons carrying the virus. The result is that each initial detection, even if isolated, should be treated as valuable data if it is going to serve as an early warning signal for in-country transmission.

The analysis of cumulative detections per site per lineage (Fig. 3) show that in most instances the time between initial and subsequent or consistent detection in Terminal 1 and 3 WW was longer compared to both the pooled aircraft sewage and municipal WW regions. This was expected when comparing municipal and terminal WW because transmission is less likely to occur within the population as the infected individuals may have been in transit. If transmission did occur within a terminal, the secondarily infected individuals may not have contributed to the sampled WW before departing. It was expected that the pooled aircraft WW would show the same lag between initial detection and subsequent detections given the transient nature of the population represented but pooled aircraft sewage did not lag to the same degree. It is possible that the passive samplers used initially were less sensitive, causing the detections that were observed after the switch to autosamplers to represent higher frequency lineage presence.

It is not typical to place confidence in single detections from WW surveillance due to the nature of the sample type, being comprised of many different individual viral particles and with each particle contributing a partially degraded genome. Frequency prediction tools are necessary to ‘detangle’ the sample into individual lineages and there are limits to the accuracy of such tools. For example, a benchmarking study which included Alcov concluded that calls at frequencies less than 5% should be interpreted with caution due to high background noise in WW samples⁴³. As a result, in municipal WW surveillance, patterns and trends in data, especially in slowly emerging lineages, lend confidence to the predictions. Initial detections in this study were predicted at 5% or higher (apart from BA.2.75 and BQ*), but there are other factors that can result in lower confidence in single-point signal spikes, such as the possibility that diagnostic mutations were missed. For this reason, in the cases (BQ*, XBB.1.5*) where the initial detection was more than 2 weeks ahead of any other detections considered in this study (other WW collected from the sites included in this analysis or Ontario clinical) the second detections should be considered as well when calculating lead time gained (bracketed values in Table 2). If WW sequencing in airports were to be used to inform pathogen impact on public health and inform policy changes in future pandemics, recommendations from data collectors to the recipients of the data would have to be well balanced to account for differences between single and multiple detections with corresponding interpretations weighted accordingly.

It is valuable as well to be able to track the potential national entry of pathogens and their emergence in the community through airport coupled with municipal WW surveillance amidst limited clinical testing to inform how detections at entry points translate to transmission in the surrounding regions³⁵. The long timespan between terminal and municipal WW detections in some instances (e.g. BQ*, Fig. 2C) suggests that in some cases variants may not have been quick to spread after initial introduction to the local population, or that they were present in airport WW due to transitioning travelers. The transmissibility and infection rates of SARS-CoV-2 differ due to many factors (i.e. by variant, age demographics, and vaccination rates of individuals), and not all viral detections directly translate to continued transmission⁵⁶. After initial detection of emerging lineages, the signal in Toronto Pearson WW generally continued to fluctuate, with higher peaks in frequency than the municipal WW from the surrounding regions. Toronto Pearson WW then would not be useful to monitor subsequent local trends in VOC dominance past the initial community introduction, which has been the key role of municipal WW testing in SARS-CoV-2 surveillance.

Unique mutations influence the ability to accurately call specific sub-lineages

XBB lineages provided an example of how frequency predictions using WW sequencing in this study were affected by availability of unique mutations, despite the ability of prediction tools including Alcov to use shared mutations as well to make lineage predictions. XBB*, as a recombinant of BA.2.10.1 and BA.2.75 sub-lineages (BJ.1 and BM.1.1.1 respectively), had a distinct mutation constellation compared to BQ* and BF* lineages, which descended from BA.5 and other BA.5 sub-lineages which were in circulation around the time of XBB* emergence in Ontario^47,57. Despite high transmissibility and associated rapid spread, XBB* was detected in Toronto Pearson WW a week ahead of the first clinical case⁴⁸. XBB.1.5 had only three nucleotide mutations (C44T, T17124C, T23018C) differentiating it from XBB.1, with only one of these resulting in an amino acid change (T23018C yields S: F486P), making it less distinct for de-mixing analyses³³. WW data collected from Toronto Pearson was unable to provide more than one early detection compared to clinical sequencing in the case of XBB.1.5. This is likely because Alcov was unable to confidently and specifically predict XBB.1.5 when it was present in samples at low abundance due to similarity to other XBB* lineages. XBB.1.5* was also prevalent in the USA before emergence in Ontario so it is likely that early transmission occurred across the Canada-USA border land crossings, which would affect lead time⁵⁸. XBB.1.16* had 12 nucleotide changes (C11750T, C11956T, T12730A, A14856G, G18703T, A22101T, C22995G, T23018C, G27915T, T28297C, A28447G, C29386T) compared to XBB.1, and XBB.1.9* had five nucleotide changes (G5720A, C12789T, T23018C, G27915T, T28297C) compared to XBB.1 which made these two groups of XBB sub-lineages slightly easier to differentiate and thus detect in WW (Table 2)³³. In addition to the differences in number of unique mutations, viral shedding across lineages/sublineages may not be the same⁵⁹. For example, one lineage might be shed for a longer period than other lineages. This in turn might affect the proportion and amplification of new variants in WW⁵⁹.

Terminal and aircraft sewage are most valuable when used together

Only Terminal 1 and 3 WW provided an average lead time prior to the PHO initial clinical detection. We would expect initial detection of emerging lineages which originate outside of Canada to be in the pooled aircraft sewage more often due to the contributing populations being more transient, but in this study emerging VOCs were more frequently detected in terminal WW first. The pooled aircraft sewage was collected using passive samplers from January 1 until October 21, 2022, much longer than the Terminals were, and many of the VOCs in this study emerged before October 2022. While passive samplers were being used, the coverage for pooled aircraft sewage samples was lower than after the auto sampler was installed (Supplementary Table ^S3), which could be contributing to the lack of lead time, though the median BOC was lower in the terminals than in the passive sampler pooled aircraft samples. The later detection and lack of average lead time of many lineages by pooled aircraft WW is also confounded by a possible reluctance of passengers to use the airplane lavatory, in favor of waiting until landing at the airport, or that the lavatories may not be needed on shorter flights⁶⁰. After October 21st 2022 when an autosampler was used for pooled aircraft sewage as well, BOC was on average much higher than in Terminals 1 and 3 (Supplementary Table ^S3) likely due to the triturator sewage being less dilute. In the examples of two XBB sub-lineages, XBB.1.9 and XBB.1.16 (Fig. 2f and d, supplementary Figures S6 and S7), which emerged in Ontario after the switch in collection method, the pooled aircraft sewage samples had earlier signals and in the case of XBB.1.16 provided the only early detection. On the other hand, in the example of the S: R346T mutation, the first detection in pooled aircraft sewage occurred within the same week as terminal first detections even though this was during the period when passive samplers were in use at the triturator. Additionally, the first detection of BF* in Terminals was in a sample collected using a passive sampler, yet it still resulted in a lead time (Supplementary Table ^S1, Fig. 2b). Whether the differences in performance between the two sample types were due to the timing of sampling method change or dilution of the sample by local staff members, both sample types will be valuable to collect in future pandemic preparedness surveillance programs as they provide insight into different population types.

Effect of sampling method, sequencing methods and sampling site

There are variations in how samples were processed in this study due to the collaboration of three sequencing laboratories to generate data for sub-sets of sampling sites. BOC was analyzed and compared in Supplementary Table ^S3. More of the locations along the SARS-CoV-2 genome where defining mutations could be found were captured in sequence data with higher BOC, which aided in detecting and distinguishing lineages. All sample groups had median BOC greater than 75% of the genome so all sites had many samples with close to complete genomes. If method, rather than sampling site, was influencing first detection location it would be expected that the University of Guelph Pooled Aircraft samples would detect VOCs first most often due to higher BOC, but this was not observed. While municipal WW did have lower average BOC all-together and first detections did not occur in municipal WW during this study, the Peel region average BOC was slightly higher than the average BOC for Terminals 1 and 3 of Toronto Pearson, which had the most frequent first detections. Samples processed by the lab at University of Waterloo provided an opportunity to analyze how the same methods yielded different results for different sampling locations as well, with average BOC differing between regions.

It is possible that the differences in sample extraction methods biased the fragments of SARS-CoV-2 viral RNA that were captured prior to ARTIC amplification and sequencing. The sequencing labs which contributed data to this study from University of Guelph and University of Waterloo used Ceres Nanotrap Microbiome A Particles to capture viral particles and other microbes from the supernatant fraction of the WW while Western University sequencing lab pelleted the sample by centrifugation. A study comparing extraction methods found that Ceres Nanotrap beads and centrifugation to extract from the solids pellet had similar human viral read yields⁶¹. The same study also found that while the centrifugation method had higher RNA yields and longer fragments than nanotrap beads, the centrifugation method did not produce complete human viral genomes from sequencing after probe-capture enrichment⁶¹. In this study the Toronto region had the lowest average BOC which is consistent with these findings though using tiled amplicon sequencing likely reduced the differences observed between extraction methods. Additionally, the Qiagen RNA extraction kit used at University of Guelph was designed for Viral RNA capture while the University of Waterloo and Western University use a more general RNA extraction kit. This could have had an influence on BOC, as samples processed through University of Guelph had higher overall BOC.

Target selection as a limitation

Pandemic preparedness through WW surveillance at airports will also require knowledge of a potentially emerging virus prior to its arrival in Canada. In the current study, tiled amplicons designed specifically to capture only SARS-CoV-2 were used to sequence the SARS-CoV-2 particles found in WW. This strategy works well for catching emerging VOCs because once the amplicons have been sequenced, the assembled genomes can be re-analyzed looking for new lineages retrospectively as new VOCs become concerning in other parts of the world. Even so, this strategy cannot identify new lineages. New lineage prediction could be accomplished using bioinformatics approaches such as non-negative matrix factorization (NMF) in the future⁶².

WW surveillance and clinical lineage reporting by GISAID or PHO reports also depend on lineage definition availability. If lineages are introduced into a community via infected individuals shortly after the initial evolution of the lineage, a lack of defined mutation constellation or lineage definition can delay identification of the lineage and therefor data availability. While clinical sequences, if near-complete, can contribute to efforts to define new lineages, WW sequencing cannot due to the complex nature of the viral community it captures.

If the surveillance is extended to other viral targets, sequencing strategies will have to be developed for each different target ahead of their potential spread into Canada for the data to be available early enough to serve as an early warning for public health decision making. This could mean developing and implementing more tiled assays for virus families that are an ongoing concern (RSV, Influenza) or those that have caused outbreaks in recent history i.e. Mpox⁶³. Alternatively, deep shotgun metagenomic sequencing could be performed to collect data without a target for detection of potentially novel viral threats. It should be noted that shotgun sequencing of viral RNA is resource-intensive as viral reads are often less abundant than bacterial reads, requiring a large amount of sequencing space and/or rigorous pre-selection or rRNA depletion to capture the viral community in detail^25,26. RT-qPCR has been used to detect a wider panel of pathogens in aircraft and airport terminal WW and could be used to detect a known, pre-determined set of pathogens in future surveillance programs, but would be limited by assay development if used to follow virus evolution²³.

The authors of this manuscript would like to thank and acknowledge the support of Public Health Agency of Canada, the Greater Toronto Airports Authority and Ontario WSI for partnering with academic labs for WW surveillance and for the support provided during the project. We would also like to thank the regional Public Health Units that provided samples for comparison and lent us their expertise during the writing of this paper; Regional Municipality of Peel, Toronto Public Health, Toronto Water, and York Region Public Health. In the Ontario WSI program, sequencing labs rely on other academic labs which perform qPCR to co-ordinate the transportation of municipal samples and ensure that all samples that need to be sequenced are forwarded in a timely fashion. As a consequence, this data could not have been collected without the attention and support provided by Dr. Mark Servos and team at the University of Waterloo, Dr. Elizabeth Edwards, Dr. Hui Peng and team at the University of Toronto and Dr. Kimberly Gilbride, Dr. Claire Oswald and team at Toronto Metropolitan University.We would like to extend our gratitude for the organization, funding and support that the Ontario Ministry of Environment Conservation and Parks (MECP) provided for this work through their leadership of the Ontario Wastewater Surveillance Initiative that has dedicated resources to the establishment and funding of a province-wide wastewater testing program. We would like to thank Bahar Aminvaziri from the MECP for her contribution to the manuscript editing process. Funding for the collection and analysis of Airport WW samples at the time of this study was provided by the Public Health Agency of Canada and Health Canada. We would like to thank Health Canada, the Public Health Agency of Canada and GTAA for the unique opportunity for research provided in their partnership with academic labs in the WSI.