The COVID-19 pandemic has highlighted the importance of genomic epidemiology in deciphering the origin and spread of SARS-CoV-2 lineages across local and global scales to aid in directing responses (Deng et al., 2020; Geoghegan et al., 2020; Gonzalez-Reiche et al., 2020; du Plessis et al., 2021; Worobey et al., 2020). Studying evolutionary relationships between SARS-CoV-2 sequences over time and space allows estimation of the relative contributions of international and domestic transmission in association with public health interventions. Thus, phylogenetics can be highly informative in evaluating the effectiveness of non-pharmaceutical interventions at curbing importations and reducing transmission.

Phylogenetic analyses are dependent upon timely generation and sharing of publicly available genetic sequences and associated metadata, such as through the Global Initiative on Sharing All Influenza Data (GISAID) platform (Khare et al., 2021; Shu and McCauley, 2017). The global availability of SARS-CoV-2 genomes has been unprecedented, such that by 15 June 2021, just over a year and a half since the first whole genome sequence was shared (de Maio, 2020; Wu et al., 2020), there were 495,159 SARS-CoV-2 sequences on GISAID representing 108,301,802 COVID-19 diagnoses (Krispin and Byrnes, 2020). Nomenclature systems to partition viral sequences sharing common mutations and recent ancestry, such as Pango lineages, have provided a dynamic system for genomic SARS-CoV-2 surveillance (O’Toole et al., 2021; Rambaut et al., 2020). Tracking the emergence, dispersal, and genomic characteristics of lineages, particularly variants of concern (VOCs) and interest (VOIs), has become critical in light of their demonstrated increased transmissibility, potential increased clinical severity, and ability to circumvent host immune responses (Lopez Bernal et al., 2021; Faria et al., 2021; Meng et al., 2021; Planas et al., 2021; Tegally et al., 2020). Therefore, identifying and characterizing predominantly circulating lineages and sublineages is a cornerstone of global epidemiological surveillance and policy.

Phylogeographic methods to infer sampled viruses’ dispersal have been widely applied to quantify SARS-CoV-2 introductions into the UK (du Plessis et al., 2021), the USA (Deng et al., 2020; Gonzalez-Reiche et al., 2020; Worobey et al., 2020; Zeller et al., 2021), Brazil (Candido et al., 2020), New Zealand (NZ) (Douglas et al., 2021; Geoghegan et al., 2020), and Europe (Hodcroft et al., 2021; Worobey et al., 2020), among others, elucidating variable epidemic dynamics. A recent review expounds further on the approaches and application of phylodynamic models towards improving our understanding of SARS-CoV-2 transmission and control (Attwood et al., 2022). In the UK, where the most sequences have been generated per case globally (Furuse, 2021), there were an estimated 1179 introductions resulting in two or more sampled descendant cases in the early epidemic (du Plessis et al., 2021). By contrast, the first wave of COVID-19 in Louisiana, USA, was primarily attributable to a single domestic introduction several weeks prior to Mardi Gras, resulting in multiple superspreader events and wide dissemination across Southern USA (Zeller et al., 2021). In NZ, where stringent border closures and lockdown measures were enacted early on, there were estimated to have been 277 introductions up to 1 July 2020, among which 19% resulted in multiple downstream cases (Geoghegan et al., 2020). Large-scale SARS-CoV-2 genomic epidemiology analyses in Canada have thus far been limited to a study on the early epidemic in the province of Quebec in which they conservatively estimated at least 500 viral introductions into Quebec by June 2020 largely attributable to the province’s spring break (Murall et al., 2021). Our analyses elaborate upon their findings at a national scale for the first and second waves of COVID-19 in Canada.

Characterization of viral importations over time can also clarify the effectiveness of public health interventions by associating inflection points in importations with drastic changes in policies such as international travel restrictions (Magalis et al., 2020). Compartmental modelling approaches have also been applied to quantify the impact of COVID-19 control measures such as social distancing, informing decisions about when to relax or increase stringency (Anderson et al., 2020; Anderson et al., 2021). Evaluating the effectiveness of these interventions is key to reacting proportionally to the risk posed by future outbreaks of SARS-CoV-2 and other zoonotic pathogens. It is unclear how changes in intervention stringency, social and travel behaviour, and circulating viral diversity affected Canadian SARS-CoV-2 transmission dynamics, particularly during a time of negligibly low immunity prior to vaccine roll-out or widespread natural infections.

The first COVID-19 case in Canada was detected on 25 January 2020 in a traveller from Wuhan to Toronto and by 5 March, the first community transmission case was identified (Press, 2021; Figure 1C). In subsequent weeks, the stringency of Canadian interventions increased rapidly, summarized by the Oxford stringency index (Hale et al., 2021). On 14 March, a travel advisory warning against all non-essential travel outside Canada was issued; on 18 March, travel restrictions on the entry of all foreign nationals (except from the USA) were implemented; 21 March travel restrictions were extended to the USA; and 24 March, a mandatory 14-day self-isolation for those returning from international travel was implemented. The stringency index during the first wave increased most rapidly on 16 March 2020, reaching the highest stringency between 1 and 3 April 2020. The first wave reached a maximum of 1108 average daily new COVID-19 diagnoses across Canada on 3 May 2020. While many of the earliest cases in Canada were attributable to A.1 and B.4 lineages in British Columbia (Figure 1—figure supplement 1), cases in March through June 2020 were primarily attributable to B.1* lineages and B.1.1* lineages in Ontario and Quebec, as well as lineage B.1.279 in Alberta in March and April. By June, in response to reductions in daily new diagnoses country-wide, restrictions slowly relaxed, resulting in a reduction of the stringency index, albeit with varying rates by province (Cameron-Blake et al., 2021). Fewer than 186 average daily new diagnoses occurred across Canada throughout summer 2020, mostly in Quebec, Ontario, Alberta, and British Columbia. Then, with the rising second wave in fall 2020, the stringency index rebounded in November 2020, although with simultaneous easing of entry exceptions for foreign nationals and quarantine shortening (Cameron-Blake et al., 2021). The second wave entered exponential growth in November 2020, as VOCs and VOIs began to displace wild-type lineages.

Figure 1 with 2 supplements see all

Download asset Open asset

The VOCs and VOIs first detected in Canada were B.1.1.7 (Alpha) on 8 November 2020, followed by P.2 on 24 November, A.23.1 on 2 December, B.1.429 on 5 December, B.1.525 (Eta) on 14 December, B.1.351 (Beta) on 19 December, B.1.427 (Epsilon) on 30 December (dates from GISAID). The winter holidays were followed by the crest of the second wave on 11 January 2021 with 3555 average daily new diagnoses in Canada, likely dampened by restrictions on travel, personal gatherings, dining, and mask use. The new year also brought first Canadian detections of B.1.526.1/.2 (Iota) on 1 January 2021, P.1 (Gamma) on 1 February, and C.37 (Lambda) on 15 February. By March 2021, the trough of the second wave gave heed to the third wave, driven almost exclusively by B.1.1.7 and P.1, the latter in Western provinces primarily. Monitoring lineage frequencies alone does not inform how many individual importations accounted for the detected cases, nor the domestic spread dynamics of individual sublineages, warranting a genomic epidemiology analysis.

A phylogenetic maximum likelihood approach was applied to estimate the timing, origins, destinations, and transmission of SARS-CoV-2 introductions in Canada from the beginning of 2020 to the end of the second wave in March 2021 for Canadian sublineages (introductions with sampled descendants) and singletons (introductions with no sampled descendants). We tested the hypothesis that international travel restrictions enacted in March 2020 effectively reduced international importations of SARS-CoV-2 into Canada, yet ongoing introductions contributed to COVID-19 persistence into early 2021, exacerbated by highly transmissible B.1.1.7 and other VOC sublineages. These analyses help to elucidate the relative contributions made by sublineages introduced before and after enactment of travel restrictions towards persistence of the Canadian SARS-CoV-2 epidemic in 2020 and early 2021, prior to the global predominance of VOCs and vaccine roll-out. Evaluating travel restrictions’ effectiveness towards reducing importations could inform the stringency and timing of future public health interventions.

Together, these analyses support the conclusion that travel restrictions and other non-pharmaceutical interventions imposed in March 2020 drastically reduced importations, preventing the expansion of subsequent sublineages. A more rapid and stringent public health response would have reduced the initial burden by preventing early sublineages from establishing widespread transmission chains that persisted throughout the first and second waves. However, wild-type (pre-VOC) sublineages introduced in summer 2020 when prevalence and immunity were low contributed the highest proportion of COVID-19 cases in the second wave, suggesting that even a low level of ongoing importations of similarly transmissible variants can contribute to viral persistence. Over time, improvements in quarantine, contact tracing, testing, and individual behavioural changes collectively contributed to smaller sublineages with shorter transmission lifespans, prior to the seeding of highly transmissible VOC sublineages. A moderate level of border porosity throughout 2020 and 2021 provided opportunities for importations of B.1.1.7 and other VOCs that outcompeted previously circulating sublineages to near extinction in the third wave.

Although the USA was highly represented in all of our subsamples as a result of its large contribution to COVID-19 cases in 2020 and high sequence availability during early months, our results suggest a greater effect than due to sampling alone. On average, the USA sequences represented 28.9% (28.7–29.2%) of total international sequences, yet accounted for 46.3% (44.0–48.7%) of all sublineages and 57.7% (55.6–59.8%) of singletons. Upon maximizing the number of Canadian sequences in the analysis, where global sequence representation was more normalized but less comprehensive, the USA sequences represented fewer of the international sequences (25.8%, 25.6–26.1%) and still accounted for 38.4% (37.0–39.8%) of sublineages and 46.4% (44.6–48.3%) of singletons. It was foreseeable that the USA would be the largest source of SARS-CoV-2 importations into Canada given its high COVID-19 prevalence throughout 2020 as well as the fact that these nations share the world’s longest land border, spanning 8890 km. Overall international arrivals into Canada declined 77.8% from 96.8 M in 2019 to 19.7 M in 2020 (Canada Statistics, 2021a). However, the number of truck drivers and crew members (air, ship, and train), only declined by 24.8%, accounting for 49% of all international arrivals after April 2020 (of whom 93% were truck drivers). Truck drivers supported the supply chain throughout COVID-19 disruptions, yet in doing so, may have inadvertently facilitated additional SARS-CoV-2 importations from the USA. International arrivals of Canadians and non-residents in 2020 were low, but not negligible, with 14.6 M and 5.1 M arrivals respectively in 2020 (Canada Statistics, 2021a). While essential workers were exempt from quarantine, all others were mandated a self-enforced 14-day quarantine upon re-entry. It was not until 22 February 2021 that a 3-day hotel quarantine was made mandatory for international arrivals by air (Government of Canada, 2021). The downstream effects of border porosity could have been better mitigated by earlier widespread availability of rapid antigen tests to inform quarantine duration (Centers for Disease Control and Prevention, 2020), and through increased resources for contact tracing, cluster response, and compliance oversight. While the majority of importations were inferred to have been into Quebec and Ontario, this could be partially attributable to differences in sequence generation and deposition. Based on provincial population sizes alone, one would expect more importations into Ontario and Quebec (populations of 14.8 M and 8.6 M) than into British Columbia (5.2 M), Alberta (4.5 M), or other provinces and territories with less than 1.5 M (Canada Statistics, 2021b).

During the first wave, early sublineages had the opportunity to become established and resulted in large transmission chains driving COVID-19 burden, consistent with findings from the UK (du Plessis et al., 2021). However, the second wave was mostly driven by cases from newly seeded sublineages. We accounted for potential confounding due to recent sublineages having been monitored for less time by stratifying sublineages by whether they were active or extinct in the previous 2 months; and by restricting the data to collection dates 3 months prior to the data download, we accounted for delay in depositing sequence data for more recent sublineages. Differences remained between early and late sublineages’ evolutionary dynamics in the presence of VOCs with heightened transmissibility and the roll-out of vaccines. The earliest introductions were typified by high sublineage lifespans (i.e. duration of transmission), which manifested as an initial increase in the mean time since importation, corroborating that outbreak control via sublineage detection and contract tracing improved over the course of the first wave. Sublineage TMRCAs typically lag importation by several days, but could precede them if multiple infections with subsequently sampled descendants occurred during travel. While accounting for sampling lag between symptoms and sample collection would get closer to the importation date, this assumes symptom onset occurred during travel, which might not always be the case. In contrast, if we assume that infected travellers who sparked domestic transmission chains were tested upon or soon after arrival, then this lag would be zero to a couple of days. If an infectious traveller was not sampled, but gave rise to an outbreak that was later sampled, then an additional lag must be considered regarding the length and number generations of unsampled viral predecessors since importation. A TMRCA could precede importation if the sublineage index case infected multiple people during travel. Where travel history is available, incorporating internal nodes or branch states for travel location would help to inform the timing as well as geographic origins of importations (Hong et al., 2021). In light of these additional sources of uncertainty, we abstained from applying an importation lag to TMRCAs, and emphasize their interpretation as the estimated time of the recent common ancestor of all sampled descendants.

Case burden in the second wave was primarily due to sublineages introduced in the summer of 2020 amid low prevalence, despite a low number of importations. This suggests that ongoing seeding events impacted viral persistence prior to the arrival of VOCs with increased transmissibility. Considering the global and temporal diversity of the SARS-CoV-2 pandemic, it is difficult to predict what variant may be introduced and subsequently disperse through a local population. While broad and longstanding restrictions against non-essential international travel is not necessarily an advisable policy in light of economic impacts, swift and stringent travel bans towards locations harbouring a high frequency of a VOC not yet identified domestically should be seriously considered to reduce the probability of seeding multiple, simultaneous outbreaks and potentially overwhelming the healthcare system. Dynamic travel bans can be detrimental from a socioeconomic perspective; therefore their implementation must be weighed relative to other non-pharmaceutical interventions to reduce transmission. As illustrated by Omicron, perhaps pervasive spread is inevitable with the most transmissible and immune evasive variants. However, reducing the number of importations sparking domestic outbreaks grants jurisdictions more time to prepare for the inevitable by ramping up vaccinations, scaling up testing and contact tracing, allocating primary care resources, and implementing non-pharmaceutical interventions. These advanced planning steps also serve to minimize the economic harm due to work absenteeism and supply chain disruption.

Ancestral trait reconstruction is sensitive to sampling bias as ancestral node states are more likely to reflect states predominant among tips (Baele et al., 2017; Hill et al., 2021; Lemey et al., 2009). By subsampling sequences proportionally to monthly relative case counts, a prior is assumed towards the probability of importation from each geography. While subsampling proportionally to cases reduces representation of geographies with the highest sequence contributions, it does not directly increase representation of geographies with sparse sequences, which may have relatively underestimated relative migration rates. Furthermore, countries and provinces differed in the extent of their contact tracing, testing, and case reporting, which adds uncertainty in the ascertainment rate underlying observed cases. Hospitalizations or COVID-related deaths could be used to impute or deconvolute the true number of people infected (Huisman et al., 2021), but this requires strong assumptions about the hospitalization and case fatality rate, both functions of the virulence of viral variants and prior immunological exposure of the population. In the present analysis, we make the simplifying assumption that case detection rates are comparable across geographies, justifying the use of case contributions to inform the subsampling process, however we recognize that additional data clarifying reason for testing could improve these estimates. The subsampling sensitivity analysis revealed how subsampling iteratively reduced overrepresentation in available provincial sequences and improves the correlation between cases and sequences. Across subsampling strategies, the relative spatiotemporal trends were robust, although exact estimates were variable, particularly for smaller values.

Incorporating individuals’ travel history or flight volume data could help to improve representation of undersampled geographies (Lemey et al., 2020; du Plessis et al., 2021). However, Bayesian phylogeography remains limited in its scalability despite improved efficiencies in the memory required (Didelot et al., 2018; du Plessis et al., 2021), and they are still limited to thousands, rather than tens or hundreds of thousand sequences. Therefore, the maximum likelihood approach that we have applied offers scalability, while taking into consideration the uncertainty of sampling bias by including multiple bootstraps and subsampling strategies. Stratified analyses of well-supported (sub)lineages in future analyses will facilitate direct comparisons of Bayesian and maximum likelihood phylogeography methodologies towards estimating importation rates.

Low sequence representation can lead to underestimates of total introductions if neither index case nor descendants were sampled, underestimates of sublineage size if not all descendants were sampled, and similarly, overestimates of the proportion of singletons, which may have been from unsampled transmission chains. Extrapolating an upper estimate of introductions is challenging in the absence of additional data. Clean genomes available in Canada prior to 1 March 2021 represented 4.2% of confirmed diagnoses (and 3.2% when 75% of Canadian sequences retained). Diagnoses were estimated to represent about 9% of total cases in Canada up to September 2020, while other geographies ranged from 5% in Italy to 99% in Qatar (Noh and Danuser, 2021). The probability of a case being detected is affected by geography (sociodemographic structure, testing capacity, and recommendations), by individual (age, contact-traced, political beliefs, co-morbidities), and by lineage (symptom severity, infectivity profile). Reason for sequencing is not always random – it could be for an outbreak investigation or to confirm VOC identity – and it varies over time by jurisdiction. As more sequences are generated and made available, we expect more descendants of previously identified sublineages than travellers or their recent contacts harbouring new sublineages or singletons. When sequencing efforts or resources are lower, travellers are a more efficient use of resources if prevalence is higher abroad than domestically, increasing the travel bias. Thus, importations do not scale linearly with sequence representation. In theory, the upper limit of importations by province could be estimated by adjusting for monthly sequence representation, case ascertainment rate, outbreak bias (ratio of probabilities of testing given infected for random versus outbreak-linked), and travel bias (ratio of probabilities of testing given infected for domestic versus travelling populations) over time, stratified by geography. More consistent inclusion of the reason for sequencing and testing in the publicly available metadata could facilitate better estimates of the extent of travel-related and outbreak-related bias. Additionally, prospective cohort studies or seroprevalence studies would ameliorate our estimate of the case ascertainment fraction.

These analyses shed light upon the natural epidemiological history of SARS-CoV-2 in the context of public health interventions and exemplify a sublineage-based method of genomic surveillance. Sharing viral genome sequences linked with the time and place of sampling in a timely manner is of utmost importance for epidemic surveillance of new and already described variants, analyses to support contact tracing, and inference of SARS-CoV-2 importation dynamics.

All the viral genetic sequences analyzed in this study were sourced from the Global initiative on sharing all influenza data (GISAID) coronavirus (CoV) database and are subject to the GISAID EpiFlu Database Access Agreement. Within this agreement, we are disallowed to distribute data to any third party other than Authorized Users as contemplated by this Agreement. In lieu, a table of accession IDs was provided in Supplementary file 1, along with an acknowledgement of all originating and submitting laboratories. Submission of this table is also a condition of use of the data under the terms and conditions of GISAID. Individuals who would like to become Authorized Users may submit an application https://www.gisaid.org/registration/register/ and agree to the Database Access Agreement. Full and subsampled alignments can be shared to Authorized Users upon request. The scripts used for data curation, inferences, analyses, and visualization are available with a reproducible example at https://github.com/AngMcL/sars-cov-2_canada_2020 (copy archived at swh:1:rev:1cb49fa43c752b1ce1003c52e43b5c47d3b6789f).

1. 1. Shu Y
   2. McCauley J

   (2017) GISAID

   GISAID: Global initiative on sharing all influenza data - coronovius (CoV).

   https://doi.org/10.2807/1560-7917.ES.2017.22.13.30494

1. 1. Anderson SC
   2. Edwards AM
   3. Yerlanov M
   4. Mulberry N
   5. Stockdale JE
   6. Iyaniwura SA
   7. Falcao RC
   8. Otterstatter MC
   9. Irvine MA
   10. Janjua NZ
   11. Coombs D
   12. Colijn C

   (2020) Quantifying the impact of covid-19 control measures using a bayesian model of physical distancing

   PLOS Computational Biology 16:e1008274.

   https://doi.org/10.1371/journal.pcbi.1008274

   • PubMed
   • Google Scholar
2. 1. Anderson SC
   2. Mulberry N
   3. Edwards AM
   4. Stockdale JE
   5. Iyaniwura SA
   6. Falcao RC
   7. Otterstatter MC
   8. Janjua NZ
   9. Coombs D
   10. Colijn C

   (2021) How much leeway is there to relax COVID-19 control measures?

   Epidemics 35:100453.

   https://doi.org/10.1016/j.epidem.2021.100453

   • PubMed
   • Google Scholar
3. 1. Attwood SW
   2. Hill SC
   3. Aanensen DM
   4. Connor TR
   5. Pybus OG

   (2022) Phylogenetic and phylodynamic approaches to understanding and combating the early SARS-CoV-2 pandemic

   Nature Reviews. Genetics pp. 1–16.

   https://doi.org/10.1038/s41576-022-00483-8

   • PubMed
   • Google Scholar
4. 1. Baele G
   2. Suchard MA
   3. Rambaut A
   4. Lemey P

   (2017) Emerging concepts of data integration in pathogen phylodynamics

   Systematic Biology 66:e47–e65.

   https://doi.org/10.1093/sysbio/syw054

   • PubMed
   • Google Scholar
5. Software

   1. Baquero OS

   (2017) ggsn: North symbols and scale bars for maps created with “ggplot2” or “ggmap.”, version 0.4.0

   R Package.

   https://CRAN.R-project.org/package=ggsn
6. Software

   1. Bedford T

   (2021) Nextstrain/ncovnextstrain/ncov, version v12

   Github.

   https://github.com/nextstrain/ncov
7. Software

   1. Bivand R
   2. Keitt T
   3. Rowlingson B

   (2017) rgdal: Bindings for the “geospatial” data abstraction library, version 1.2-16

   R Package.

   https://CRAN.R-project.org/package=rgdal
8. Software

   1. Bivand RS
   2. Lewin-Koh N

   (2018a) Tools for handling spatial objects, version 0.9-4

   R package.

   https://CRAN.R-project.org/package=maptools
9. Software

   1. Bivand R
   2. Rundel C

   (2018b) rgeos: Interface to geometry engine - open source ('GEOS’), version 0.4-1

   R Package.

   https://CRAN.R-project.org/package=rgeos
10. 1. Brunson J

    (2020) Ggalluvial: Layered grammar for alluvial plots

    Journal of Open Source Software 5:e2017.

    https://doi.org/10.21105/joss.02017

    • Google Scholar
11. Website

    1. Cameron-Blake E
    2. Breton C
    3. Sim P
    4. Tatlow H
    5. Hale T
    6. Wood A
    7. Smith J
    8. Sawatsky J
    9. Parsons Z
    10. Tyson K

    (2021) Variation in the Canadian provincial and territorial responses to COVID-19

    Accessed July 11, 2022.

    https://www.bsg.ox.ac.uk/covidtracker
12. Website

    1. Canada Statistics

    (2019) 2016 census - boundary files

    Statcan.gc.ca. Accessed November 13, 2019.

    https://www12.statcan.gc.ca/census-recensement/2011/geo/bound-limit/bound-limit-2016-eng.cfm
13. Website

    1. Canada Statistics

    (2021a) Table 24-10-0041-01 International travellers entering or returning to Canada, by type of transport

    150.statcan.gc.ca. Accessed July 14, 2022.

    https://doi.org/10.25318/2410004101-eng
14. Website

    1. Canada Statistics

    (2021b) Canada’s population clock (real-time model)

    150.statcan.gc.ca. Accessed June 22, 2022.

    https://www150.statcan.gc.ca/n1/pub/71-607-x/71-607-x2018005-eng.htm
15. Website

    1. Canadian Institute for Health Information

    (2021) Covid-19 intervention timeline in canada - data tables

    Canadian Institute for Health Information. Accessed June 9, 2022.

    https://www.cihi.ca/en/covid-19-intervention-timeline-in-canada
16. 1. Candido DS
    2. Claro IM
    3. de Jesus JG
    4. Souza WM
    5. Moreira FRR
    6. Dellicour S
    7. Mellan TA
    8. du Plessis L
    9. Pereira RHM
    10. Sales FCS
    11. Manuli ER
    12. Thézé J
    13. Almeida L
    14. Menezes MT
    15. Voloch CM
    16. Fumagalli MJ
    17. Coletti TM
    18. da Silva CAM
    19. Ramundo MS
    20. Amorim MR
    21. Hoeltgebaum HH
    22. Mishra S
    23. Gill MS
    24. Carvalho LM
    25. Buss LF
    26. Prete CA
    27. Ashworth J
    28. Nakaya HI
    29. Peixoto PS
    30. Brady OJ
    31. Nicholls SM
    32. Tanuri A
    33. Rossi ÁD
    34. Braga CKV
    35. Gerber AL
    36. de C Guimarães AP
    37. Gaburo N
    38. Alencar CS
    39. Ferreira ACS
    40. Lima CX
    41. Levi JE
    42. Granato C
    43. Ferreira GM
    44. Francisco RS
    45. Granja F
    46. Garcia MT
    47. Moretti ML
    48. Perroud MW
    49. Castiñeiras T
    50. Lazari CS
    51. Hill SC
    52. de Souza Santos AA
    53. Simeoni CL
    54. Forato J
    55. Sposito AC
    56. Schreiber AZ
    57. Santos MNN
    58. de Sá CZ
    59. Souza RP
    60. Resende-Moreira LC
    61. Teixeira MM
    62. Hubner J
    63. Leme PAF
    64. Moreira RG
    65. Nogueira ML
    66. Ferguson NM
    67. Costa SF
    68. Proenca-Modena JL
    69. Vasconcelos ATR
    70. Bhatt S
    71. Lemey P
    72. Wu CH
    73. Rambaut A
    74. Loman NJ
    75. Aguiar RS
    76. Pybus OG
    77. Sabino EC
    78. Faria NR
    79. Brazil-UK Centre for Arbovirus Discovery, Diagnosis, Genomics and Epidemiology (CADDE) Genomic Network

    (2020) Evolution and epidemic spread of SARS-CoV-2 in Brazil

    Science 369:1255–1260.

    https://doi.org/10.1126/science.abd2161

    • PubMed
    • Google Scholar
17. Website

    1. Centers for Disease Control and Prevention

    (2020) Science brief- options to reduce quarantine of persons with sars-cov-2 infection using s monitoring and diagnostic testing

    Pagefreezer. Accessed May 11, 2022.

    https://www.cdc.gov/coronavirus/2019-ncov/more/scientific-brief-options-to-reduce-quarantine.html
18. Software

    1. Couch DR

    (2021) Broom: Convert statistical objects into tidy tibbles r package version 0.7.6, version 0.7.6

    R Package.

    https://CRAN.R-project.org/package=broom
19. Website

    1. de Maio N

    (2020) Novel 2019 coronavirus genome

    Virological.org. Accessed February 28, 2020.

    https://virological.org/t/novel-2019-coronavirus-genome/319
20. 1. Deng X
    2. Gu W
    3. Federman S
    4. du Plessis L
    5. Pybus OG
    6. Faria NR
    7. Wang C
    8. Yu G
    9. Bushnell B
    10. Pan CY
    11. Guevara H
    12. Sotomayor-Gonzalez A
    13. Zorn K
    14. Gopez A
    15. Servellita V
    16. Hsu E
    17. Miller S
    18. Bedford T
    19. Greninger AL
    20. Roychoudhury P
    21. Starita LM
    22. Famulare M
    23. Chu HY
    24. Shendure J
    25. Jerome KR
    26. Anderson C
    27. Gangavarapu K
    28. Zeller M
    29. Spencer E
    30. Andersen KG
    31. MacCannell D
    32. Paden CR
    33. Li Y
    34. Zhang J
    35. Tong S
    36. Armstrong G
    37. Morrow S
    38. Willis M
    39. Matyas BT
    40. Mase S
    41. Kasirye O
    42. Park M
    43. Masinde G
    44. Chan C
    45. Yu AT
    46. Chai SJ
    47. Villarino E
    48. Bonin B
    49. Wadford DA
    50. Chiu CY

    (2020) Genomic surveillance reveals multiple introductions of SARS-CoV-2 into Northern California

    Science 369:582–587.

    https://doi.org/10.1126/science.abb9263

    • PubMed
    • Google Scholar
21. 1. Didelot X
    2. Croucher NJ
    3. Bentley SD
    4. Harris SR
    5. Wilson DJ

    (2018) Bayesian inference of ancestral dates on bacterial phylogenetic trees

    Nucleic Acids Research 46:e134.

    https://doi.org/10.1093/nar/gky783

    • PubMed
    • Google Scholar
22. Software

    1. Dinno A

    (2017) Dunn.Test: Dunn’s test of multiple comparisons using rank sums, version 1.3.5

    R Package.

    https://CRAN.R-project.org/package=dunn.test
23. 1. Douglas J
    2. Geoghegan JL
    3. Hadfield J
    4. Bouckaert R
    5. Storey M
    6. Ren X
    7. de Ligt J
    8. French N
    9. Welch D

    (2021) Real-time genomics for tracking severe acute respiratory syndrome coronavirus 2 border incursions after virus elimination, new zealand

    Emerging Infectious Diseases 27:2361–2368.

    https://doi.org/10.3201/eid2709.211097

    • PubMed
    • Google Scholar
24. 1. du Plessis L
    2. McCrone JT
    3. Zarebski AE
    4. Hill V
    5. Ruis C
    6. Gutierrez B
    7. Raghwani J
    8. Ashworth J
    9. Colquhoun R
    10. Connor TR
    11. Faria NR
    12. Jackson B
    13. Loman NJ
    14. O’Toole Á
    15. Nicholls SM
    16. Parag KV
    17. Scher E
    18. Vasylyeva TI
    19. Volz EM
    20. Watts A
    21. Bogoch II
    22. Khan K
    23. Aanensen DM
    24. Kraemer MUG
    25. Rambaut A
    26. Pybus OG
    27. COVID-19 Genomics UK (COG-UK) Consortium

    (2021) Establishment and lineage dynamics of the SARS-CoV-2 epidemic in the UK

    Science 371:708–712.

    https://doi.org/10.1126/science.abf2946

    • PubMed
    • Google Scholar
25. 1. Faria NR
    2. Mellan TA
    3. Whittaker C
    4. Claro IM
    5. Candido D
    6. Mishra S
    7. Crispim MAE
    8. Sales FCS
    9. Hawryluk I
    10. McCrone JT
    11. Hulswit RJG
    12. Franco LAM
    13. Ramundo MS
    14. de Jesus JG
    15. Andrade PS
    16. Coletti TM
    17. Ferreira GM
    18. Silva CAM
    19. Manuli ER
    20. Pereira RHM
    21. Peixoto PS
    22. Kraemer MUG
    23. Gaburo N
    24. Camilo C
    25. Hoeltgebaum H
    26. Souza WM
    27. Rocha EC
    28. de Souza LM
    29. de Pinho MC
    30. Araujo LJT
    31. Malta FSV
    32. de Lima AB
    33. Silva J
    34. Zauli DAG
    35. Ferreira A
    36. Schnekenberg RP
    37. Laydon DJ
    38. Walker PGT
    39. Schlüter HM
    40. Dos Santos ALP
    41. Vidal MS
    42. Del Caro VS
    43. Filho RMF
    44. Dos Santos HM
    45. Aguiar RS
    46. Proença-Modena JL
    47. Nelson B
    48. Hay JA
    49. Monod M
    50. Miscouridou X
    51. Coupland H
    52. Sonabend R
    53. Vollmer M
    54. Gandy A
    55. Prete CA
    56. Nascimento VH
    57. Suchard MA
    58. Bowden TA
    59. Pond SLK
    60. Wu CH
    61. Ratmann O
    62. Ferguson NM
    63. Dye C
    64. Loman NJ
    65. Lemey P
    66. Rambaut A
    67. Fraiji NA
    68. Carvalho M
    69. Pybus OG
    70. Flaxman S
    71. Bhatt S
    72. Sabino EC

    (2021) Genomics and epidemiology of the P.1 SARS-CoV-2 lineage in Manaus, Brazil

    Science 372:815–821.

    https://doi.org/10.1126/science.abh2644

    • PubMed
    • Google Scholar
26. 1. Furuse Y

    (2021) Genomic sequencing effort for SARS-CoV-2 by country during the pandemic

    International Journal of Infectious Diseases 103:305–307.

    https://doi.org/10.1016/j.ijid.2020.12.034

    • PubMed
    • Google Scholar
27. 1. Geoghegan JL
    2. Ren X
    3. Storey M
    4. Hadfield J
    5. Jelley L
    6. Jefferies S
    7. Sherwood J
    8. Paine S
    9. Huang S
    10. Douglas J
    11. Mendes FK
    12. Sporle A
    13. Baker MG
    14. Murdoch DR
    15. French N
    16. Simpson CR
    17. Welch D
    18. Drummond AJ
    19. Holmes EC
    20. Duchêne S
    21. de Ligt J

    (2020) Genomic epidemiology reveals transmission patterns and dynamics of SARS-CoV-2 in Aotearoa New Zealand

    Nature Communications 11:1–7.

    https://doi.org/10.1038/s41467-020-20235-8

    • PubMed
    • Google Scholar
28. 1. Gonzalez-Reiche AS
    2. Hernandez MM
    3. Sullivan MJ
    4. Ciferri B
    5. Alshammary H
    6. Obla A
    7. Fabre S
    8. Kleiner G
    9. Polanco J
    10. Khan Z
    11. Alburquerque B
    12. van de Guchte A
    13. Dutta J
    14. Francoeur N
    15. Melo BS
    16. Oussenko I
    17. Deikus G
    18. Soto J
    19. Sridhar SH
    20. Wang YC
    21. Twyman K
    22. Kasarskis A
    23. Altman DR
    24. Smith M
    25. Sebra R
    26. Aberg J
    27. Krammer F
    28. García-Sastre A
    29. Luksza M
    30. Patel G
    31. Paniz-Mondolfi A
    32. Gitman M
    33. Sordillo EM
    34. Simon V
    35. van Bakel H

    (2020) Introductions and early spread of SARS-CoV-2 in the New York City area

    Science 369:297–301.

    https://doi.org/10.1126/science.abc1917

    • PubMed
    • Google Scholar
29. Website

    1. Government of Canada

    (2021) Government of canada expands restrictions to international travel by land and air.Pdf

    Government of Canada. Accessed February 22, 2021.

    https://www.canada.ca/en/public-health/news/2021/02/government-of-canada-expands-restrictions-to-international-travel-by-land-and-air.html
30. 1. Grolemund G
    2. Wickham H

    (2011) Dates and times made easy with lubridate

    Journal of Statistical Software 40:ev040i03.

    https://doi.org/10.18637/jss.v040.i03

    • Google Scholar
31. 1. Hadfield J
    2. Megill C
    3. Bell SM
    4. Huddleston J
    5. Potter B
    6. Callender C
    7. Sagulenko P
    8. Bedford T
    9. Neher RA

    (2018) Nextstrain: Real-time tracking of pathogen evolutionNextstrain: real-time tracking of pathogen evolution

    Bioinformatics 34:4121–4123.

    https://doi.org/10.1093/bioinformatics/bty407

    • PubMed
    • Google Scholar
32. 1. Hale T
    2. Angrist N
    3. Goldszmidt R
    4. Kira B
    5. Petherick A
    6. Phillips T
    7. Webster S
    8. Cameron-Blake E
    9. Hallas L
    10. Majumdar S
    11. Tatlow H

    (2021) A global panel database of pandemic policies (Oxford COVID-19 Government Response Tracker)

    Nature Human Behaviour 5:529–538.

    https://doi.org/10.1038/s41562-021-01079-8

    • PubMed
    • Google Scholar
33. Software

    1. Henry L
    2. Wickham H
    3. Chang W

    (2020) gGstance: Horizontal “ggplot2” components, version 0.3.5

    R Package.

    https://CRAN.R-project.org/package=ggstance
34. 1. Hill V
    2. Ruis C
    3. Bajaj S
    4. Pybus OG
    5. Kraemer MUG

    (2021) Progress and challenges in virus genomic epidemiology

    Trends in Parasitology 37:1038–1049.

    https://doi.org/10.1016/j.pt.2021.08.007

    • PubMed
    • Google Scholar
35. 1. Hodcroft EB
    2. Zuber M
    3. Nadeau S
    4. Vaughan TG
    5. Crawford KHD
    6. Althaus CL
    7. Reichmuth ML
    8. Bowen JE
    9. Walls AC
    10. Corti D
    11. Bloom JD
    12. Veesler D
    13. Mateo D
    14. Hernando A
    15. Comas I
    16. González-Candelas F
    17. Stadler T
    18. Neher RA
    19. SeqCOVID-SPAIN consortium

    (2021) Spread of a SARS-CoV-2 variant through Europe in the summer of 2020

    Nature 595:707–712.

    https://doi.org/10.1038/s41586-021-03677-y

    • PubMed
    • Google Scholar
36. 1. Hong SL
    2. Lemey P
    3. Suchard MA
    4. Baele G

    (2021) Bayesian phylogeographic analysis incorporating predictors and individual travel histories in beast

    Current Protocols 1:e98.

    https://doi.org/10.1002/cpz1.98

    • PubMed
    • Google Scholar
37. Preprint

    1. Huisman JS
    2. Scire J
    3. Angst DC
    4. Li J
    5. Neher RA
    6. Maathuis MH
    7. Bonhoeffer S
    8. Stadler T

    (2021) Estimation and Worldwide Monitoring of the Effective Reproductive Number of SARS-Cov-2Estimation and Worldwide Monitoring of the Effective Reproductive Number of SARS-CoV-2

    medRxiv.

    https://doi.org/10.1101/2020.11.26.20239368

    • Google Scholar
38. 1. Kalia K
    2. Saberwal G
    3. Sharma G

    (2021) The lag in SARS-CoV-2 genome submissions to GISAID

    Nature Biotechnology 39:1058–1060.

    https://doi.org/10.1038/s41587-021-01040-0

    • PubMed
    • Google Scholar
39. 1. Khare S
    2. Gurry C
    3. Freitas L
    4. Schultz MB
    5. Bach G
    6. Diallo A
    7. Akite N
    8. Ho J
    9. Lee RT
    10. Yeo W
    11. Curation Team GC
    12. Maurer-Stroh S

    (2021) Gisaid’s role in pandemic response

    China CDC Weekly 3:1049–1051.

    https://doi.org/10.46234/ccdcw2021.255

    • PubMed
    • Google Scholar
40. Software

    1. Krispin R
    2. Byrnes J

    (2020) Coronavirus: The 2019 novel coronavirus COVID-19 (2019-ncov) datasetcoronavirus: The 2019 Novel Coronavirus COVID-19 (2019-nCoV) Dataset, version 0.3.0

    Github.

    https://github.com/RamiKrispin/coronavirus
41. Software

    1. Lanfear R

    (2020) A global phylogeny of SARS-CoV-2 sequences from GISAID, version 3.0

    Github.

    https://github.com/roblanf/sarscov2phylo
42. 1. Lemey P
    2. Rambaut A
    3. Drummond AJ
    4. Suchard MA

    (2009) Bayesian phylogeography finds its roots

    PLOS Computational Biology 5:e1000520.

    https://doi.org/10.1371/journal.pcbi.1000520

    • PubMed
    • Google Scholar
43. 1. Lemey P
    2. Hong SL
    3. Hill V
    4. Baele G
    5. Poletto C
    6. Colizza V
    7. O’Toole Á
    8. McCrone JT
    9. Andersen KG
    10. Worobey M
    11. Nelson MI
    12. Rambaut A
    13. Suchard MA

    (2020) Accommodating individual travel history and unsampled diversity in bayesian phylogeographic inference of sars-cov-2

    Nature Communications 11:1–14.

    https://doi.org/10.1038/s41467-020-18877-9

    • PubMed
    • Google Scholar
44. 1. Li H

    (2018) Minimap2: Pairwise alignment for nucleotide sequencesMinimap2: pairwise alignment for nucleotide sequences

    Bioinformatics 34:3094–3100.

    https://doi.org/10.1093/bioinformatics/bty191

    • PubMed
    • Google Scholar
45. 1. Lopez Bernal J
    2. Andrews N
    3. Gower C
    4. Gallagher E
    5. Simmons R
    6. Thelwall S
    7. Stowe J
    8. Tessier E
    9. Groves N
    10. Dabrera G
    11. Myers R
    12. Campbell CNJ
    13. Amirthalingam G
    14. Edmunds M
    15. Zambon M
    16. Brown KE
    17. Hopkins S
    18. Chand M
    19. Ramsay M

    (2021) Effectiveness of covid-19 vaccines against the b.1.617.2 (delta) variant

    The New England Journal of Medicine 385:585–594.

    https://doi.org/10.1056/NEJMoa2108891

    • PubMed
    • Google Scholar
46. 1. Lu J
    2. du Plessis L
    3. Liu Z
    4. Hill V
    5. Kang M
    6. Lin H
    7. Sun J
    8. François S
    9. Kraemer MUG
    10. Faria NR
    11. McCrone JT
    12. Peng J
    13. Xiong Q
    14. Yuan R
    15. Zeng L
    16. Zhou P
    17. Liang C
    18. Yi L
    19. Liu J
    20. Xiao J
    21. Hu J
    22. Liu T
    23. Ma W
    24. Li W
    25. Su J
    26. Zheng H
    27. Peng B
    28. Fang S
    29. Su W
    30. Li K
    31. Sun R
    32. Bai R
    33. Tang X
    34. Liang M
    35. Quick J
    36. Song T
    37. Rambaut A
    38. Loman N
    39. Raghwani J
    40. Pybus OG
    41. Ke C

    (2020) Genomic epidemiology of sars-cov-2 in guangdong province, china

    Cell 181:997–1003.

    https://doi.org/10.1016/j.cell.2020.04.023

    • PubMed
    • Google Scholar
47. Preprint

    1. Magalis BR
    2. Ramirez-Mata A
    3. Zhukova A
    4. Mavian C
    5. Marini S
    6. Lemoine F
    7. Prosperi M
    8. Gascuel O
    9. Salemi M

    (2020) Differing Impacts of Global and Regional Responses on SARS-Cov-2 Transmission Cluster DynamicsDiffering Impacts of Global and Regional Responses on SARS-CoV-2 Transmission Cluster Dynamics

    bioRxiv.

    https://doi.org/10.1101/2020.11.06.370999

    • Google Scholar
48. Preprint

    1. McLaughlin A
    2. Montoya V
    3. Miller RL
    4. Mordecai GJ
    5. Worobey M
    6. Poon AFY
    7. Joy JB

    (2021) Early and Ongoing Importations of SARS-Cov-2 in CanadaEarly and Ongoing Importations of SARS-CoV-2 in Canada

    medRxiv.

    https://doi.org/10.1101/2021.04.09.21255131

    • Google Scholar
49. 1. Meng B
    2. Kemp SA
    3. Papa G
    4. Datir R
    5. Ferreira I
    6. Marelli S
    7. Harvey WT
    8. Lytras S
    9. Mohamed A
    10. Gallo G
    11. Thakur N
    12. Collier DA
    13. Mlcochova P
    14. Duncan LM
    15. Carabelli AM
    16. Kenyon JC
    17. Lever AM
    18. De Marco A
    19. Saliba C
    20. Culap K
    21. Cameroni E
    22. Matheson NJ
    23. Piccoli L
    24. Corti D
    25. James LC
    26. Robertson DL
    27. Bailey D
    28. Gupta RK
    29. COVID-19 Genomics UK (COG-UK) Consortium

    (2021) Recurrent emergence of sars-cov-2 spike deletion h69/v70 and its role in the alpha variant b.1.1.7

    Cell Reports 35:109292.

    https://doi.org/10.1016/j.celrep.2021.109292

    • PubMed
    • Google Scholar
50. 1. Minh BQ
    2. Schmidt HA
    3. Chernomor O
    4. Schrempf D
    5. Woodhams MD
    6. von Haeseler A
    7. Lanfear R

    (2020) Iq-tree 2: New models and efficient methods for phylogenetic inference in the genomic era

    Molecular Biology and Evolution 37:1530–1534.

    https://doi.org/10.1093/molbev/msaa015

    • PubMed
    • Google Scholar
51. 1. Moshiri N
    2. Robinson P

    (2021) Viralmsa: Massively scalable reference-guided multiple sequence alignment of viral genomesViralMSA: massively scalable reference-guided multiple sequence alignment of viral genomes

    Bioinformatics 37:714–716.

    https://doi.org/10.1093/bioinformatics/btaa743

    • Google Scholar
52. 1. Murall CL
    2. Fournier E
    3. Galvez JH
    4. N’Guessan A
    5. Reiling SJ
    6. Quirion PO
    7. Naderi S
    8. Roy AM
    9. Chen SH
    10. Stretenowich P
    11. Bourgey M
    12. Bujold D
    13. Gregoire R
    14. Lepage P
    15. St-Cyr J
    16. Willet P
    17. Dion R
    18. Charest H
    19. Lathrop M
    20. Roger M
    21. Bourque G
    22. Ragoussis J
    23. Shapiro BJ
    24. Moreira S

    (2021) A small number of early introductions seeded widespread transmission of SARS-CoV-2 in Québec, Canada

    Genome Medicine 13:169.

    https://doi.org/10.1186/s13073-021-00986-9

    • PubMed
    • Google Scholar
53. Software

    1. Neuwirth E

    (2014) Rcolorbrewer: Colorbrewer palettes, version 1.1-2

    R Package.

    https://CRAN.R-project.org/package=RColorBrewer
54. 1. Noh J
    2. Danuser G

    (2021) Estimation of the fraction of COVID-19 infected people in U.S. states and countries worldwide

    PLOS ONE 16:e0246772.

    https://doi.org/10.1371/journal.pone.0246772

    • PubMed
    • Google Scholar
55. 1. O’Toole Á
    2. Scher E
    3. Underwood A
    4. Jackson B
    5. Hill V
    6. McCrone JT
    7. Colquhoun R
    8. Ruis C
    9. Abu-Dahab K
    10. Taylor B
    11. Yeats C
    12. du Plessis L
    13. Maloney D
    14. Medd N
    15. Attwood SW
    16. Aanensen DM
    17. Holmes EC
    18. Pybus OG
    19. Rambaut A

    (2021) Assignment of epidemiological lineages in an emerging pandemic using the pangolin tool

    Virus Evolution 7:eveab064.

    https://doi.org/10.1093/ve/veab064

    • PubMed
    • Google Scholar
56. Software

    1. Pagès H
    2. Aboyoun P
    3. Gentleman R
    4. DebRoy S

    (2020) Biostrings: Efficient manipulation of biological strings, version 2.56.0

    R Package.

    https://bioconductor.org/packages/release/bioc/html/Biostrings.html
57. 1. Paradis E
    2. Claude J
    3. Strimmer K

    (2004) Ape: Analyses of phylogenetics and evolution in R languageAPE: Analyses of Phylogenetics and Evolution in R language

    Bioinformatics 20:289–290.

    https://doi.org/10.1093/bioinformatics/btg412

    • PubMed
    • Google Scholar
58. 1. Paradis E
    2. Schliep K

    (2019) Ape 5.0: An environment for modern phylogenetics and evolutionary analyses in rape 5.0: an environment for modern phylogenetics and evolutionary analyses in R

    Bioinformatics 35:526–528.

    https://doi.org/10.1093/bioinformatics/bty633

    • PubMed
    • Google Scholar
59. 1. Planas D
    2. Veyer D
    3. Baidaliuk A
    4. Staropoli I
    5. Guivel-Benhassine F
    6. Rajah MM
    7. Planchais C
    8. Porrot F
    9. Robillard N
    10. Puech J
    11. Prot M
    12. Gallais F
    13. Gantner P
    14. Velay A
    15. Le Guen J
    16. Kassis-Chikhani N
    17. Edriss D
    18. Belec L
    19. Seve A
    20. Courtellemont L
    21. Péré H
    22. Hocqueloux L
    23. Fafi-Kremer S
    24. Prazuck T
    25. Mouquet H
    26. Bruel T
    27. Simon-Lorière E
    28. Rey FA
    29. Schwartz O

    (2021) Reduced sensitivity of SARS-CoV-2 variant Delta to antibody neutralization

    Nature 596:276–280.

    https://doi.org/10.1038/s41586-021-03777-9

    • PubMed
    • Google Scholar
60. Website

    1. Press C

    (2021) A timeline of COVID-19 in Canada

    Nationalpost.com. Accessed January 24, 2021.

    https://nationalpost.com/pmn/news-pmn/canada-news-pmn/a-timeline-of-covid-19-in-canada
61. 1. Price MN
    2. Dehal PS
    3. Arkin AP

    (2010) FastTree 2--approximately maximum-likelihood trees for large alignments

    PLOS ONE 5:e9490.

    https://doi.org/10.1371/journal.pone.0009490

    • PubMed
    • Google Scholar
62. Website

    1. Public Health Agency of Canada

    (2021) Coronavirus disease 2019 (COVID-19): Epidemiology update

    Government of Canada. Accessed April 1, 2021.

    https://health-infobase.canada.ca/covid-19/epidemiological-summary-covid-19-cases.html?stat=num&measure=active#a2
63. 1. Rambaut A
    2. Lam TT
    3. Max Carvalho L
    4. Pybus OG

    (2016) Exploring the temporal structure of heterochronous sequences using tempest (formerly path-o-gen)

    Virus Evolution 2:evew007.

    https://doi.org/10.1093/ve/vew007

    • PubMed
    • Google Scholar
64. Website

    1. Rambaut A

    (2020) Phylodynamic analysis | 176 genomes |

    Virological.org. Accessed May 6, 2020.

    https://virological.org/t/phylodynamic-analysis-176-genomes-6-mar-2020/356/1
65. 1. Rambaut A
    2. Holmes EC
    3. O’Toole Á
    4. Hill V
    5. McCrone JT
    6. Ruis C
    7. du Plessis L
    8. Pybus OG

    (2020) A dynamic nomenclature proposal for SARS-CoV-2 lineages to assist genomic epidemiology

    Nature Microbiology 5:1403–1407.

    https://doi.org/10.1038/s41564-020-0770-5

    • PubMed
    • Google Scholar
66. 1. Revell LJ

    (2012) Phytools: An r package for phylogenetic comparative biology (and other things)

    Methods in Ecology and Evolution 3:217–223.

    https://doi.org/10.1111/j.2041-210X.2011.00169.x

    • Google Scholar
67. 1. Schliep K
    2. Potts AJ
    3. Morrison DA
    4. Grimm GW
    5. Fitzjohn R

    (2017) Intertwining phylogenetic trees and networks

    Methods in Ecology and Evolution 8:1212–1220.

    https://doi.org/10.1111/2041-210X.12760

    • Google Scholar
68. 1. Shu Y
    2. McCauley J

    (2017) GISAID: Global initiative on sharing all influenza data - from vision to reality

    Euro Surveillance 22:30494.

    https://doi.org/10.2807/1560-7917.ES.2017.22.13.30494

    • PubMed
    • Google Scholar
69. Software

    1. Slowikowski K

    (2021) Ggrepel: Automatically position non-overlapping text labels with “ggplot2.”, version 0.9.1

    R package.

    https://CRAN.R-project.org/package=ggrepel
70. Preprint

    1. Tegally H
    2. Wilkinson E
    3. Giovanetti M
    4. Iranzadeh A
    5. Fonseca V
    6. Giandhari J
    7. Doolabh D
    8. Pillay S
    9. San EJ
    10. Msomi N
    11. Mlisana K
    12. von Gottberg A
    13. Walaza S
    14. Allam M
    15. Ismail A
    16. Mohale T
    17. Glass AJ
    18. Engelbrecht S
    19. Van Zyl G
    20. Preiser W
    21. Petruccione F
    22. Sigal A
    23. Hardie D
    24. Marais G
    25. Hsiao M
    26. Korsman S
    27. Davies MA
    28. Tyers L
    29. Mudau I
    30. York D
    31. Maslo C
    32. Goedhals D
    33. Abrahams S
    34. Laguda-Akingba O
    35. Alisoltani-Dehkordi A
    36. Godzik A
    37. Wibmer CK
    38. Sewell BT
    39. Lourenço J
    40. Alcantara LCJ
    41. Pond SLK
    42. Weaver S
    43. Martin D
    44. Lessells RJ
    45. Bhiman JN
    46. Williamson C
    47. de Oliveira T

    (2020) Emergence and Rapid Spread of a New Severe Acute Respiratory Syndrome-Related Coronavirus 2 (SARS-Cov-2) Lineage with Multiple Spike Mutations in South AfricaEmergence and Rapid Spread of a New Severe Acute Respiratory Syndrome-Related Coronavirus 2 (SARS-CoV-2) Lineage with Multiple Spike Mutations in South Africa

    medRxiv.

    https://doi.org/10.1101/2020.12.21.20248640

    • Google Scholar
71. 1. To TH
    2. Jung M
    3. Lycett S
    4. Gascuel O

    (2016) Fast dating using least-squares criteria and algorithms

    Systematic Biology 65:82–97.

    https://doi.org/10.1093/sysbio/syv068

    • PubMed
    • Google Scholar
72. Website

    1. United States Census Bureau

    (2019) Cartographic boundary files

    Map. Accessed April 7, 2022.

    https://www.census.gov/geographies/mapping-files/time-series/geo/cartographic-boundary.html
73. Book

    1. Venables WN
    2. Ripley BD

    (2002) Modern Applied Statistics with S

    Springer.

    https://doi.org/10.1007/978-0-387-21706-2

    • Google Scholar
74. 1. Wickham H
    2. Averick M
    3. Bryan J
    4. Chang W
    5. McGowan L
    6. François R
    7. Grolemund G
    8. Hayes A
    9. Henry L
    10. Hester J
    11. Kuhn M
    12. Pedersen T
    13. Miller E
    14. Bache S
    15. Müller K
    16. Ooms J
    17. Robinson D
    18. Seidel D
    19. Spinu V
    20. Takahashi K
    21. Vaughan D
    22. Wilke C
    23. Woo K
    24. Yutani H

    (2019) Welcome to the tidyverseWelcome to the Tidyverse

    Journal of Open Source Software 4:e1686.

    https://doi.org/10.21105/joss.01686

    • Google Scholar
75. Software

    1. Wilke CO

    (2020) Cowplot: Streamlined plot theme and plot annotations for “ggplot2.”, version 1.1.1

    R package.

    https://CRAN.R-project.org/package=cowplot
76. Software

    1. Wilke CO

    (2021) Ggridges: Ridgeline plots in “ggplot2.”, version 0.5.3

    R package.

    https://CRAN.R-project.org/package=ggridges
77. Software

    1. Wilkins D

    (2021) Treemapify: Draw treemaps in “ggplot2.”, version 2.5.5

    R package.

    https://CRAN.R-project.org/package=treemapify
78. 1. Worobey M
    2. Pekar J
    3. Larsen BB
    4. Nelson MI
    5. Hill V
    6. Joy JB
    7. Rambaut A
    8. Suchard MA
    9. Wertheim JO
    10. Lemey P

    (2020) The emergence of SARS-CoV-2 in Europe and North America

    Science 370:564–570.

    https://doi.org/10.1126/science.abc8169

    • PubMed
    • Google Scholar
79. 1. Wu F
    2. Zhao S
    3. Yu B
    4. Chen YM
    5. Wang W
    6. Song ZG
    7. Hu Y
    8. Tao ZW
    9. Tian JH
    10. Pei YY
    11. Yuan ML
    12. Zhang YL
    13. Dai FH
    14. Liu Y
    15. Wang QM
    16. Zheng JJ
    17. Xu L
    18. Holmes EC
    19. Zhang YZ

    (2020) A new coronavirus associated with human respiratory disease in China

    Nature 579:265–269.

    https://doi.org/10.1038/s41586-020-2008-3

    • PubMed
    • Google Scholar
80. 1. Yu G
    2. Smith DK
    3. Zhu H
    4. Guan Y
    5. Lam TTY

    (2016) Ggtree: An rpackage for visualization and annotation of phylogenetic trees with their covariates and other associated data

    Methods in Ecology and Evolution 8:28–36.

    https://doi.org/10.1111/2041-210x.12628

    • Google Scholar
81. Software

    1. Yu G

    (2020) Ggplotify: Convert plot to “grob” or “ggplot” object, version 0.0.5

    R package.

    https://CRAN.R-project.org/package=ggplotify
82. 1. Zeileis A
    2. Grothendieck G

    (2005) Zoo: S3 infrastructure for regular and irregular time series

    Journal of Statistical Software 14:e06.

    https://doi.org/10.18637/jss.v014.i06

    • Google Scholar
83. 1. Zeller M
    2. Gangavarapu K
    3. Anderson C
    4. Smither AR
    5. Vanchiere JA
    6. Rose R
    7. Snyder DJ
    8. Dudas G
    9. Watts A
    10. Matteson NL
    11. Robles-Sikisaka R
    12. Marshall M
    13. Feehan AK
    14. Sabino-Santos G
    15. Bell-Kareem AR
    16. Hughes LD
    17. Alkuzweny M
    18. Snarski P
    19. Garcia-Diaz J
    20. Scott RS
    21. Melnik LI
    22. Klitting R
    23. McGraw M
    24. Belda-Ferre P
    25. DeHoff P
    26. Sathe S
    27. Marotz C
    28. Grubaugh ND
    29. Nolan DJ
    30. Drouin AC
    31. Genemaras KJ
    32. Chao K
    33. Topol S
    34. Spencer E
    35. Nicholson L
    36. Aigner S
    37. Yeo GW
    38. Farnaes L
    39. Hobbs CA
    40. Laurent LC
    41. Knight R
    42. Hodcroft EB
    43. Khan K
    44. Fusco DN
    45. Cooper VS
    46. Lemey P
    47. Gardner L
    48. Lamers SL
    49. Kamil JP
    50. Garry RF
    51. Suchard MA
    52. Andersen KG

    (2021) Emergence of an early SARS-CoV-2 epidemic in the United States

    Cell 184:4939–4952.

    https://doi.org/10.1016/j.cell.2021.07.030

    • PubMed
    • Google Scholar

Decision letter after peer review:

Thank you for submitting your article "Early introductions of SARS-CoV-2 sublineages into Canada drove the 2020 epidemic" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by David Serwadda as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Bernardo Gutierrez (Reviewer #3).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

All reviewers thought this study represented a potentially valuable contribution. They also all expressed concern about the impact of potentially biased sampling on the conclusions. We recognize that this is a general challenge in the field, and not all papers in it have necessarily faced this level of scrutiny. Because the conclusions of this paper are quantitative, it seems important to quantify the uncertainty (in this case, to model misspecification/sampling bias) accurately.

The reviewers outline nuances in their individual commentaries below, but there are effectively three general issues:

1) It seems particularly important to understand (e.g., as a sensitivity analysis) the impact of not downsampling the samples from Canada.

2) Reviewers also asked whether the case-based scaling of sequences is appropriate. This relates to a more general problem raised by the highly uneven sampling by country and region. Downsampling so that no region is overrepresented with respect to its relative infection load would leave zero samples for analysis, since so many areas are missing sequences entirely. Perhaps again performing some kind of sensitivity analysis or referencing a more thorough investigation of this problem, how does uneven representation impact estimated migration rates?

3) The reviewers were also concerned with the estimate of the total number of introductions.

Effectively, all the reviews are requesting a more careful treatment of potential problems with data collection (it is not random!), both in the analyses and in the text.

Reviewer #1 (Recommendations for the authors):

A few other statistical questions are raised in the recommendations for the authors below.

Statistics/methodology

1. p. 5, l. 3: Clearer logic in the "Subsampling" section would help. What is the justification for subsampling proportional to case number, which is obviously quite biased by surveillance effort, perhaps even more biased than the number of sequences? Why not subsample proportional to countries' population size or excess mortality, adjusted for age? It would be helpful to describe too how much subsampling occurred and that it was necessary for computational reasons. What kinds of biases are still present? For instance, it would be good to point out somewhere the regions/periods that are so poorly sampled (e.g., in many lower income countries) that the inference could be biased to favor endemic transmission, assuming my intuition is correct.

2. p. 8, ll. 6-8: Could you explain the logic of how the total number of introductions was inferred given the estimated rate at which cases were sequenced (and cases detected)? Should the estimated number of introductions scale linearly?

3. p. 12, ll. 6-9 (Figure 4C): Does this model properly adjust for the typical delay between sampling members of a sublineage? Sublineages with more recent tMRCA should have fewer descendants simply because less time has passed. It is not clear the stratification takes care of this. There might also be a decrease in the number of descendants with time if the total number of sublineages is increasing and surveillance does not scale with sublineage richness. The text reads "To elucidate the relative contributions of early and late sublineages,….," but it is unclear exactly what contribution is being measured here.

4. p. 10, ll. 10-11 (Figure 4D): I think we could see an increase in days since importation over time even if more recently imported lineages were generally more common ("dominant"). The maximum possible age of the oldest lineage is increasing with time. The sentence says this "suggests" that early introductions "increasingly dominated" the epidemic. Maybe it would be easier to demonstrate this through specific examples or a measure of persistence time and a proxy of prevalence.

5. p. 18, ll. 6-8: "Sequences from individuals with travel history to Iran or Italy before June 2020 were recategorized as having been sampled in the country of travel." Why was this not done for individuals with travel history from other countries?

Interpretation

1. Abstract: "Rapid implementation of stringent border controls and quarantine could have diminished the Canadian COVID-19 burden by curtailing the spread of early introductions." I question how much stringent border controls, unless defined very severely as prohibiting the introduction of, e.g., more than ten lineages, could have diminished the COVID-19 burden. Once any lineage is established and transmitting in the community, additional introductions have negligible impact on prevalence or the "final size" of an epidemic. Quarantine, by reducing the effective reproductive number, can reduce burden. I worry that lack of precision on the impact of introductions on burden can lead to the sort of knee-jerk travel bans we see with Omicron.

2. p. 3, ll. 7-8: "can illuminate"? Deciphering the rates of imported and domestic transmission is not guaranteed.

3. p. 8, l. 8: A 42% infection detection rate (as cases) struck me as rather high. The cited reference appears to be for something else.

4. p. 10, ll. 1-3: What is the purpose of mentioning the number of global descendants when it is known to be such an underestimate?

5. Figure 4. It would be useful to reference 4A and 4B specifically in the text.

6. p. 15, ll. 4-6: How can we separate the impacts of behavioral changes (NPIs) from seasonality on sublineage size and duration?

7. p. 15, ll. 22-24: Where is the mandatory 14-day quarantine coming from? Given the incubation period and duration of shedding, I know many who think quarantine should be shorter (if testing were not feasible). With sensitive testing, there is no reason to keep people past a week if they are consistently negative. I think it's important not to imply the traditional durations used for quarantine are especially grounded or optimized.

Potential enhancements

It would be interesting to see how the importation rates from different countries correlate with air travel to/from those countries. I realize this is difficult given connecting flights, but perhaps the data exist somewhere.

p. 16, ll. 11-14: Ideally these nonstationary sampling probabilities by province (and really for all countries) would be built directly into the likelihood.

Reviewer #2 (Recommendations for the authors):

1. Line 17-21 argue genomic epidemiology is useful to track and characterize VOI/VOC, but the sequences analyzed in this manuscript are largely pre-VOC so I think it makes sense to re-focus and expand this discussion to include commentary on the importance of evaluating the success of public health interventions, travel restrictions, etc.

2. I think outlining a hypothesis at the end of the intro is useful, but as it is currently written the hypothesis is very long and difficult to understand. Can the authors restate and clarify the hypothesis into a single sentence?

3. As the authors discuss, their results and those of all genomic epidemiological studies are very sensitive to sequence quality, availability, distribution, and sampling paradigms. They reduce the impact of unequal global sequencing efforts by subsampling global sequences proportionally to their monthly case counts, excluding Canadian sequences. I understand the rationale here, but worry that inclusion of all available Canadian sequences (which are not proportional to case count trends as is clear in Figure 1D and 1E) skewed their results. Perhaps one way to assess the robustness of the results would be to repeat the main analyses using a subsampled (proportion to case counts) Canadian dataset and compare results? Relatedly, even case counts could be skewed depending on testing availability and accuracy -- could the authors include a brief commentary on the availability, accessibility, and capacity of testing centers in Canada during the first year of the pandemic?

- it is also clear to see in Figure 1C and 1F that US sequences are overrepresented in this dataset compared to case counts -- do the authors think this is inflating the proportion of US imports? (I agree it is clear most imports came from the US)

- the Alberta and Maritimes outbreaks in spring of 2020 appear very striking (large outbreaks with lots of person-to-person spread and very few descendants) in Figure. However, this interpretation is fraught when you consider Figure 1E and 1D because there are lots of available sequences during this time compared to the total case count vs fall of 2020 where there is actually a higher case burden and very few available sequences. I think many readers will get hung up on this so maybe the authors could address it briefly in the text?

- While the color scheme is lovely, it is difficult to distinguish between some provinces particularly in Figure 2

4. Lines 6-10. I think it is a big leap to estimate the total possible introductions based on sequences representative of 1.1% of total cases.

5. It might be very informative to take a more careful look at rates of viral spread (R0) over time in addition to patterns of viral importation. This could help contextualize the effect of public health interventions on domestic spread. The authors bring up quarantine, contact tracing, testing, and individual behavior change in the discussion -- do the authors think the rate of international imports is the best measure to assess these interventions? I think R0 would be just as informative!

6. I am skeptical the results on reductions in sublineage size, increasing tMRC, and perhaps even sources of domestic transmission will persist if the Canadian sequences are subsampled proportionally to case counts. Can the authors please assess the robustness of these results against a subsampled dataset?

7. I agree with and appreciate the authors' discussion of the importance data sharing, including a limited metadata set to aid in the rapid interpretation and use of sequence data for important public health decisions.

8. Although this is a stylistic preference, I believe a cleaned up version of Figure 7 with the Canada's key public health interventions and COVID-19 outbreaks along with a map of Canada's provinces (already a supplemental) and case counts during the study period would be very useful to include as an early, primary figure.

Reviewer #3 (Recommendations for the authors):

My main concern lies on the explicit naming of sources of importation following the application of an epidemiologically-based approach to downsample the publicly available genome sequences from GISAID. While it is mentioned that ten subsamples were generated, it is not clear to me that the ancestral state reconstruction was evaluated on these ten subsamples and how the results compare between them. The uncertainty of the ASR itself was incorporated for the fixed phylogeny through an ML approach (ape::ace()), but I am missing the analysis of the robustness of the results across different subsamples.

Furthermore, even if these analyses where performed, it is unclear whether a different subsampling approach could lead to the identification of a different profile of source countries, and if so, which subsampling approach would be more accurate. The accuracy of the identified countries of origin (and therefore accuracy of the proposed downsampling approach) would be validated if a complementary data source showed similar patterns: an expected importation index (EII ~ Incidence * no. of travellers to Canada), self-reported travel histories of patients with confirmed COVID-19, etc. If these complementary data sources are unavailable, I would recommend reducing the emphasis on the identification of countries of origin (or constraining it to regions where you're confident that the monthly reported cases are reliable – the USA would be an obvious example of this), or adding a clear discussion on these limitations. I would argue that this goes beyond the already mentioned idea that downsampling does not allow to increase the numbers of sequences from some unsampled times/locations (page 16, lines 5-20), and that it can also reduce the probability of phylogenetically identifying true importations from countries with low case prevalence.

On other topics: I am unsure about the estimation of total number of introductions presented in page 8 (lines 6-10). Do you have any particular reason why you would expect that importations scale linearly with sampling/representation? Given the overall dynamics, wouldn't it be more probable that higher sampling intensities would add unsampled sequences to the larger sublineages in the country (proportional to the size of the lineages) as well as uncovering unsampled importations? The rationale behind these calculations requires some more detailed explanations.

Regarding the domestic transmission, I was wondering if you had access to finer geographical resolution beyond province-level? Provinces are particularly large geographical areas, and providing some insights (rather than a full, high-resolution phylogeographic re-analysis of all the data) into the sampling within provinces could add to the manuscript. This is particularly the case for the analysis of international importations to Quebec and Ontario: where these sublineages first detected in large urban areas or not necessarily? Is the sampling within provinces representative of the cases reported in urban and rural areas? Are there any differential dynamics between larger urban areas and more remote locations?

Finally, it is quite interesting to see that the peak of TMRCAs nicely coincide with the maximum stringency on March 21 (Figure 4A,B) – even more so given that the TMRCA of a sublineage is more likely representing something akin to the first transmission event within said sublineage rather than the importation event per se. It would be interesting if the TMRCAs could be contextualised to show their distribution after accounting for sampling lag (i.e., the time between symptom onset and sample collection). This could be estimated from sequence metadata or approximated from non-genomic epidemiological data, but it would give a better idea of the times during which the infections that lead to these sublineages started.

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Genomic epidemiology of the first two waves of SARS-CoV-2 in Canada" for further consideration by eLife. Your revised article has been evaluated by Sara Sawyer (Senior Editor) and a Reviewing Editor.

The manuscript has been dramatically improved---the sensitivity analyses are immensely helpful---but there are some remaining issues that need to be addressed, as outlined below. The reviewing editor is concerned that some of the policy suggestions do not follow the evidence presented.

1. Abstract: "… suggesting travel restrictions and quarantine must be sustained to fully curtail COVID-19 burden." I am guessing most readers will likely interpret "to fully curtail COVID-19 burden" as to get the burden near zero. Is this what the authors intended? It affects the logic. The results here show that travel restrictions were usually associated with at least several importations per week. Mathematically, we know that the burden of COVID-19 does not really depend on the number of importations if R ~ 1 or more---all you need is at least one endemic lineage to maintain the prevalence/burden dictated by R. Maybe the authors are referring to conditions when R<1 and travel restrictions are even greater than what they have been historically, but this is unclear. It is also unclear precisely what the authors mean by "quarantine." Usually, it refers to the isolation of contacts of potential cases (not the cases themselves), and sometimes it refers more generally to lockdowns. Neither policy seems particularly worth calling out here. Why not say "sustained interventions to lower transmission (R<1)"? It seems unjustified for this study to single out quarantine (either def) as the necessary complementary intervention when so many other approaches have been shown more effective. I'm dwelling on the wording because these kinds of conclusions, especially in an abstract, can have an outsized impact on policy.

2. Abstract: "…restrictions that reduce the probability of importations are most effective during periods of low domestic prevalence and low or waning immunity." Waning immunity means R is rising, which, all things equal reduces extinction probabilities and increases the success of invading lineages (e.g., Otto and Whitlock 1997). Without more precision about R and how effectiveness is being defined, it's not clear restrictions should be more impactful under these conditions. I do see the positive correlation between "periods of waning immunity" and possible ban-assisted extinction, but only because immunity is most often waning when there has recently been an epidemic, and the recently boosted immune protection is relatively high. I worry again about the implications of these claims, which are not exactly supported by the research here.

3. Relatedly, the main text describes evaluating the relative contribution of (reducing) importations to (reducing) *burden* (e.g., ll. 190-192, ll. 605-607, 691-701), but this analysis is not really performed here. Such analysis would require careful estimation of R and prevalence and a good enough transmission model to evaluate extinction probabilities and outbreak sizes (in the case of R<1) over time, and it would have to account for clonal interference. Put differently, the fact that imported lineages caused x% of cases does not mean that the case count would've been (100-x)% without those importations. The fact that imported lineages displaced resident lineages does not mean the resident lineages could not have caused comparable infections had the imported lineages been stopped at the border, assuming the imported lineages had no fitness differences. Imperfect travel bans that allow some lineages to invade can at best slightly (at measured rates) change the timing of an epidemic (ll. 704-707), but usually not its size or severity. The authors write, "Although blind travel bans may not be beneficial from a socioeconomic perspective, ultimately governments need to protect citizens and dynamic travel bans are one of the few tools available." This is surprising considering the number of broader tools at our disposal (e.g., rapid testing, vaccination, ventilation/filtration) that can slow importations AND endemic spread, as well as the observed ineffectiveness of travel bans at stopping spread outside of islands and China. The authors IMO have not demonstrated how much bans could delay importation and allow time to "plan", or if investments to reduce transmission generally are more effective. I suggest these claims be carefully reconsidered.

I apologize for the wordiness. I would be happy to bolster these points (maybe with math) if they are unclear. I am concerned that too much is being extrapolated about the possible and actual epidemiological impacts of travel bans from the measured rates of importation.

You may wish to consider some suggestions to improve readability from another reviewer:

4. In line 424, I'm not entirely sure what the authors mean by "the first wave of singletons had a higher amplitude than the second wave". By 'amplitude', do they refer to the range of dates when singletons were identified between both waves? Is this relative to a feature of the importations of sublineages (given that this sentence talks about this comparison between sublineage and singleton dynamics)? Some clarification would be helpful.

5. The Discussion addresses a lot of the reviewer comments thoroughly and systematically, which was very helpful, but could be streamlined for the readers. I think that, for example, the two paragraphs that go from line 711 to line 741 could be reduced to a single paragraph discussing 1. sampling bias in phylogeography, 2. the specific effects of overrepresentation of some locations and how they addressed it with subsampling and the sensitivity analyses, 3. effects of upsampled lineages and under-sampled locations and how this can theoretically be overcome, 4. the broader challenges of scaling up Bayesian phylogenetics (this last point could be a separate paragraph, but unrelated to sampling bias).

Essential revisions:

All reviewers thought this study represented a potentially valuable contribution. They also all expressed concern about the impact of potentially biased sampling on the conclusions. We recognize that this is a general challenge in the field, and not all papers in it have necessarily faced this level of scrutiny. Because the conclusions of this paper are quantitative, it seems important to quantify the uncertainty (in this case, to model misspecification/sampling bias) accurately.

We are very grateful to the reviewers and editors for their time in reviewing our manuscript and for the many helpful comments, suggestions and clarifying questions, which we have incorporated in the revised manuscript. Most notably, we have assessed the impact of sample bias on importation rate estimates though an empirical subsampling sensitivity analysis. Additionally, we have updated the entire analysis with a more recently downloaded dataset, which included about four times the number of Canadian sequences, and an expanded study period up until March 2021 to include the entirety of the second wave. To accommodate more Canadian sequences, we have extended the global subsampling strategy to subsample Canadian sequences proportionally to provincial case counts. To address the issues of potential sampling bias on relative importations into provinces over time, we explored multiple subsampling strategies in which 25%, 50%, 75%, or 100% of Canadian sequences were retained, with the remainder of sequences from the global set up to 50,000 total. We found that under all the subsampling strategies, our conclusions were robust with regards to the impact of the travel restrictions, the temporal dynamics of importation rates, and the relative contributions of geographies towards sublineages in the first and second waves.

The reviewers outline nuances in their individual commentaries below, but there are effectively three general issues:

1) It seems particularly important to understand (e.g., as a sensitivity analysis) the impact of not downsampling the samples from Canada.

We thank the reviewers for identifying the importance of this limitation. We have addressed the impact of subsampling the Canadian sequences relative to global context sequences in a sensitivity analysis.

Methods

“Subsampling sensitivity analysis

To evaluate the effect of differing levels of Canadian sequence inclusion relative to global representation, a subsampling sensitivity analysis was conducted comparing four subsampling strategies with 25% (n = 9,184), 50% (n = 18,368), 75% (n = 27,552), or 100% (n = 36,736) of available clean Canadian sequences subsampled proportionally to monthly provincial new diagnoses and the remainder of global sequences subsampled proportionally to monthly global new diagnoses up to 50,000 total (Figure 6). […] The sensitivity of the sublineage characterization models including detection lag over time, number of descendants over time stratified by active, and sublineage lifespan stratified by active were also compared across strategies.”

Results

“Subsampling sensitivity analysis

To render phylogenetic inferences computationally feasible, reduce overrepresentation of geographies with more sequences per case, and evaluate uncertainty in the effect of sampling on importation rates, sequences were subsampled proportionally to either provinces or countries’ contributions to monthly case counts in Canada and globally for multiple bootstraps. […] Across all strategies, there were significant reductions of sublineage size and lifespan over time.”

2) Reviewers also asked whether the case-based scaling of sequences is appropriate. This relates to a more general problem raised by the highly uneven sampling by country and region. Downsampling so that no region is overrepresented with respect to its relative infection load would leave zero samples for analysis, since so many areas are missing sequences entirely. Perhaps again performing some kind of sensitivity analysis or referencing a more thorough investigation of this problem, how does uneven representation impact estimated migration rates?

We thank the reviewers for acknowledging this key issue in interpreting phylogeographic migration rates. We believe this issue has been partially addressed with the empirical subsampling sensitivity analysis. Sampling bias has been further elaborated upon in the Discussion, with additional citations on the sensitivity of discrete trait analysis to sampling bias.

“Ancestral trait reconstruction is sensitive to sampling bias as ancestral node states are more likely to reflect states predominant among tips (Baele et al., 2016; Hill et al., 2021; Lemey et al., 2009). […] Stratified analyses of well-supported (sub)lineages in future analyses will facilitate direct comparisons of Bayesian and maximum likelihood phylogeography methods employed here.”

3) The reviewers were also concerned with the estimate of the total number of introductions.

We gratefully acknowledge the reviewers’ concern regarding the rough extrapolation of an upper limit on the total number of introductions. Rather than estimate the upper limit, we have elaborated upon the type of data required to attempt this with any shred of confidence.

“Low sequence representation can lead to underestimates of total introductions if neither index case or descendants were sampled, underestimates of sublineage size if not all descendants were sampled, and similarly, overestimates of the proportion of singletons, which may have been from unsampled transmission chains. […] Additionally, prospective cohort studies or seroprevalence studies would ameliorate our estimate of the case ascertainment fraction.”

Effectively, all the reviews are requesting a more careful treatment of potential problems with data collection (it is not random!), both in the analyses and in the text.

Reviewer #1 (Recommendations for the authors):

A few other statistical questions are raised in the recommendations for the authors below.

Statistics/methodology

1. p. 5, l. 3: Clearer logic in the "Subsampling" section would help. What is the justification for subsampling proportional to case number, which is obviously quite biased by surveillance effort, perhaps even more biased than the number of sequences? Why not subsample proportional to countries' population size or excess mortality, adjusted for age?

We acknowledge the reviewer’s concern regarding the effect of differences in surveillance upon case proportional subsampling, however we maintain our perspective that the extent of bias is greater among available sequences than it is among detected cases. The supplementary figures comparing sequences per case demonstrate the utility of this subsampling method in improving the correlation of cases to sequences. While it would be possible to sample proportionally to population size, this exerts a prior belief that infected cases should be equally distributed globally and that a larger population necessarily results in a larger risk at any given time. An example corroborating this would be China, which contributed most cases in January and February 2020, but very few relative cases in late 2020 and early 2021, despite a large population size. Alternatively, COVID deaths or excess mortality could have been used, but are subject to biases of their own. Definitions and identification of COVID deaths differed over time and space, and while excess mortality would account for missed deaths, it would also include excess deaths due to other indirect causes (delayed surgeries due to COVID, for instance). We have added to the Discussion elaborating upon the justification of this method and the sensitivity analyses we conducted regarding the assumptions.

“Ancestral trait reconstruction is sensitive to sampling bias as ancestral node states are more likely to reflect states predominant among tips (Baele et al., 2016; Hill et al., 2021; Lemey et al., 2009). […] The subsampling sensitivity analysis revealed how subsampling iteratively reduced overrepresentation in available provincial sequences and improves the correlation between cases and sequences.”

It would be helpful to describe too how much subsampling occurred and that it was necessary for computational reasons. What kinds of biases are still present? For instance, it would be good to point out somewhere the regions/periods that are so poorly sampled (e.g., in many lower income countries) that the inference could be biased to favor endemic transmission, assuming my intuition is correct.

We have clarified the extent of subsampling by geography by sampling strategy.

The general effect of lingering sampling biases was elaborated upon in the Discussion (see excerpt above) and more specifically in regards to the overrepresentation of the USA:

“Although the USA was highly represented in all of our subsamples as a result of its large contribution to COVID-19 cases in 2020 and high sequence availability during early months, our results suggest a greater effect than due to sampling alone. On average, the USA sequences represented 28.9% (28.7 – 29.2%) of total international sequences, yet accounted for 46.3% (44.0 – 48.7%) of all sublineages and 57.7% (55.6 – 59.8%) of singletons. Upon maximizing the number of Canadian sequences in the analysis, where global sequence representation was more normalized but less comprehensive, the USA sequences represented fewer of the international sequences (25.8%, 25.6 – 26.1%) and still accounted for 38.4% (37.0 – 39.8%) of sublineages and 46.4% (44.6 – 48.3%) of singletons.”

The reviewer posited that the inference could be biased to favor endemic transmission, by which we assume they meant the inference would be more likely to detect local or domestic transmission sources than international sources of sampled Canadian viral diversity. This is possible in instances where there was an absence of sequences sharing a recent common ancestor with Canadian sequences from the true sublineage origin location, which itself was seeded by a previous ancestor with descendants in Canada. In the subsampling sensitivity analysis, where Canadian sequences were differentially represented in the data, the proportion of transmission attributable to international sources was comparable, although with fewer total estimated importations when too few Canadian sequences were included (25% strategy) or too few global sequences were included (100% strategy). In the 75% strategy, there were nearly the same number of Canadian and global context sequences, which if we assume a well-mixed population, confers a similar probability of detecting domestic versus international transmission sources. In reality, the global sample is not randomly drawn from a well-mixed viral pool, despite our best efforts. Naturally, we will miss detecting an importation from a geography with no sequences available. We comment upon the potential of advances in methodology by Phillipe Lemey and colleagues towards incorporating unsampled nodes in phylogeographic inference to mitigating this issue.

“Incorporating individuals’ travel history or flight volume data could help to improve representation of undersampled geographies (Lemey et al., 2020; du Plessis et al., 2021). However, Bayesian phylogeography remains limited in its scalability despite improved efficiencies in the memory required (Didelot et al., 2018; du Plessis et al., 2021), and they are still limited to thousands, rather than tens or hundreds of thousand sequences. Therefore, the maximum likelihood approach that we have applied offers scalability, while taking into consideration the uncertainty of sampling bias by including ten bootstraps for four subsampling strategies. Random sampling permits an exploration of the effect of country and Canadian province representation over time on the results. In comparing strategies, the relative spatiotemporal trends were robust, while raw estimates were more variable, as were smaller rate estimates. Stratified analyses of well-supported (sub)lineages in future analyses will facilitate direct comparisons of Bayesian and maximum likelihood phylogeography methods employed here.”

2. p. 8, ll. 6-8: Could you explain the logic of how the total number of introductions was inferred given the estimated rate at which cases were sequenced (and cases detected)? Should the estimated number of introductions scale linearly?

We gratefully acknowledge the reviewer’s concern regarding the rough extrapolation of an upper limit on the total number of introductions. Rather than estimate the upper limit, we have elaborated upon the type of data required to attempt this with any shred of confidence.

“Low sequence representation can lead to underestimates of total introductions if neither index case or descendants were sampled, underestimates of sublineage size if not all descendants were sampled, and similarly, overestimates of the proportion of singletons, which may have been from unsampled transmission chains. […] Additionally, prospective cohort studies or seroprevalence studies would ameliorate our estimate of the case ascertainment fraction.”

3. p. 12, ll. 6-9 (Figure 4C): Does this model properly adjust for the typical delay between sampling members of a sublineage? Sublineages with more recent tMRCA should have fewer descendants simply because less time has passed. It is not clear the stratification takes care of this. There might also be a decrease in the number of descendants with time if the total number of sublineages is increasing and surveillance does not scale with sublineage richness. The text reads "To elucidate the relative contributions of early and late sublineages,….," but it is unclear exactly what contribution is being measured here.

In the updated analysis, we included data up to March 2021 that had been deposited by June 2021. This permitted us to include the entire first and second waves, reducing the effect of differential delays in observing descendants. While the stratification of sublineages on the basis of having any sampled descendants in the previous two months does not entirely adjust for the delay in deposition, it improves our confidence in asserting trends among extinct sublineages. In regards to whether we might observe a decrease in the number of descendants with time if the total number of sublineages is increasing but surveillance does not scale with sublineage richness, this assumes that the testing capacity is surpassed. Until the testing capacity in a given region is reached and surpassed, the proportion of a sublineage’s descendants that are tested and sequenced should remain relatively constant, notwithstanding differences in contact tracing implementation. Although testing capacity was exceeded later in 2021, this did not affect data from the first two waves as much.

“During the first wave, early sublineages had the opportunity to become established and resulted in large transmission chains driving COVID-19 burden, consistent with findings from the UK (du Plessis et al., 2021). However, the second wave was mostly driven by cases from newly seeded sublineages. We accounted for potential confounding due to recent sublineages having been monitored for less time by stratifying sublineages by whether they were active or extinct in the previous two months; and by restricting the data to collection dates three months prior to the data download, we accounted for delay in depositing sequence data for more recent sublineages. Differences remained between early and late sublineages’ evolutionary dynamics in the presence of VOCs with heightened transmissibility and the roll-out of vaccines.”

4. p. 10, ll. 10-11 (Figure 4D): I think we could see an increase in days since importation over time even if more recently imported lineages were generally more common ("dominant"). The maximum possible age of the oldest lineage is increasing with time. The sentence says this "suggests" that early introductions "increasingly dominated" the epidemic. Maybe it would be easier to demonstrate this through specific examples or a measure of persistence time and a proxy of prevalence.

We agree with the reviewer that this requires clarification and validation. Similar to persistence time or duration of transmission, we have modelled sublineage lifespan, the time period between the TMRCA and last sample date, and stratified sublineages activity based on whether there were any samples in the last two months or if they appear to have gone extinct. By updating the analysis with data available up to 15 June 2021, but only investigating data sampled prior to 1 March, the effect of deposition delay on right truncation was reduced.

5. p. 18, ll. 6-8: "Sequences from individuals with travel history to Iran or Italy before June 2020 were recategorized as having been sampled in the country of travel." Why was this not done for individuals with travel history from other countries?

There were very few Canadian sequences with travel histories available (n=10). Besides sequences associated with Iran or Italy, these included one of each with travel history from Wuhan, Hong Kong, and an international cruise ship. Although we could have re-assigned the state of Canadian sequences with travel history to Wuhan and Hong Kong, the majority of sequences available during January and February 2020 were from China anyways. Furthermore, re-assigning tip state prevents us from identifying these sequences as introductions, so it is only a practice we would recommend for the most severely undersampled countries, which Iran and Italy were during the early pandemic.

Interpretation

1. Abstract: "Rapid implementation of stringent border controls and quarantine could have diminished the Canadian COVID-19 burden by curtailing the spread of early introductions." I question how much stringent border controls, unless defined very severely as prohibiting the introduction of, e.g., more than ten lineages, could have diminished the COVID-19 burden. Once any lineage is established and transmitting in the community, additional introductions have negligible impact on prevalence or the "final size" of an epidemic. Quarantine, by reducing the effective reproductive number, can reduce burden. I worry that lack of precision on the impact of introductions on burden can lead to the sort of knee-jerk travel bans we see with Omicron.

Thank you for your perspective on this subject, which remains an active discussion. We have elaborated upon this dialogue in the Discussion.

“The burden of cases in the second wave was primarily due to sublineages introduced in the summer of 2020, despite a low number of importations and case burden at that time. […] These advanced planning steps also serve to minimize the economic harm due to work absenteeism and supply chain disruption.”

2. p. 3, ll. 7-8: "can illuminate"? Deciphering the rates of imported and domestic transmission is not guaranteed.

Thank you for this suggestion. We have changed the wording.

“Studying evolutionary relationships between SARS-CoV-2 sequences over time and space allows estimation of the relative contributions of international and domestic transmission in association with public health interventions.”

3. p. 8, l. 8: A 42% infection detection rate (as cases) struck me as rather high. The cited reference appears to be for something else.

Thank you for this point. We have identified a more appropriate and robust reference for the case ascertainment rate in Canada in 2020.

“Clean genomes available in Canada prior to 1 March 2021 represented 4.2% of confirmed diagnoses (and 3.2% when 75% of Canadian sequences retained). Diagnoses were estimated to represent about 9% of total cases in Canada up to September 2020, while other geographies ranged from 5% in Italy to 99% in Qatar (Noh & Danuser, 2021).”

4. p. 10, ll. 1-3: What is the purpose of mentioning the number of global descendants when it is known to be such an underestimate?

Although it will naturally be an underestimate as a consequence of subsampling and incomplete case ascertainment, we have included it as a way of showing that migration flows both ways and that Canadian sublineages have been exported out of the country. While absolute numbers cannot be extracted, the relative representation of countries among sampled descendants can be interpreted with some uncertainty.

5. Figure 4. It would be useful to reference 4A and 4B specifically in the text.

Please note that the figure numbering has changed during our revisions. What were previously figure 4A and B are now Figure 5A and Figure 5—figure supplement 1. They are referenced in the text here:

“Following the entry restriction for foreign nationals in mid-March, there was a reduction, but not an elimination, of the proportion of transmission events attributable to the USA or other international sources across all provinces (Figure 5A). Prior to the restrictions, provinces had a mean of 26% of transmission events with any international origins, ranging from 17% (14-20%) in Saskatchewan to 33% (31-35%) in British Columbia. In April 2020, following travel restrictions, all provinces had decreased proportions of international transmission sources, yet Ontario and Quebec retained relatively high mean number of transmission events attributable to international sources, at 115 (105-125) and 95 (82-109) (Figure 5—figure supplement 1). This may suggest slower implementation or compliance with quarantine guidelines in these provinces.”

6. p. 15, ll. 4-6: How can we separate the impacts of behavioral changes (NPIs) from seasonality on sublineage size and duration?

We agree with the reviewer that seasonality is an important confounding factor to consider in associating differences in transmission with the implementation of NPIs or behavioral changes –factors such as temperature, humidity, and socializing inside in poorly ventilated areas can all affect transmission, and thus sublineage size and duration. In this 2020 analysis, it was not possible to disentangle these effects over the span of one year. If might be tempting to compare to 2021, but the emergence and domination of VOCs make them hard to compare. Within 2020, one might consider comparing the effect of a given NPI, such as restrictions on international travel, on sublineage size and duration in different seasons in Northern and Southern hemisphere countries. However, the challenge here is that countries different in the timing and strength of their NPIs and had different dominant sublineages circulating at any given time. Although the relationship with sublineage size and duration with time may be confounded by seasonality, the reduction in importations during Spring 2020 was too rapid and drastic to have occurred as a result of differences in season.

7. p. 15, ll. 22-24: Where is the mandatory 14-day quarantine coming from? Given the incubation period and duration of shedding, I know many who think quarantine should be shorter (if testing were not feasible). With sensitive testing, there is no reason to keep people past a week if they are consistently negative. I think it's important not to imply the traditional durations used for quarantine are especially grounded or optimized.

We have reworded the discussion to not necessarily recommend a 14-day quarantine.

“The downstream effects of border porosity could have been better mitigated by earlier widespread availability of rapid antigen tests to inform quarantine duration (CDC, 2020), and through increased resources for contact tracing, cluster response, and compliance oversight.”

Potential enhancements

It would be interesting to see how the importation rates from different countries correlate with air travel to/from those countries. I realize this is difficult given connecting flights, but perhaps the data exist somewhere.

We agree with the reviewer that this would be an interesting and challenging comparison of expected versus observed importation rates. This is similar to what was undertaken by du Plessis et al., where both prevalence and flight volume were considered towards an estimated Importation Intensity. If we had access to volumes of arrival flights or number of passengers into Canada by country for 2020, we could make this comparison. If there were more importations than expected, this may suggest systematic differences in quarantine or testing oversight based on where people arrived from. However, flight data does not incorporate land crossings, which we expect accounted for a large number of importations from the USA. We discussed the potential role played by travel across the land border towards importations from the USA.

“It was foreseeable that the USA would be the largest source of SARS-CoV-2 importations into Canada given its high COVID-19 prevalence throughout 2020 as well as the fact that these nations share the world’s longest land border, spanning 8,890 km. Overall international arrivals into Canada declined 77.8% from 96.8M in 2019 to 19.7M in 2020 (Statistics Canada, 2021a). However, the number of truck drivers and other crew members only declined by 24.8%, accounting for 49% of all international arrivals after April 2020 (of whom 93% were truck drivers). Truck drivers supported the supply chain throughout COVID-19 disruptions, yet in doing so, may have inadvertently facilitated additional SARS-CoV-2 importations from the USA. International arrivals of Canadians and non-residents in 2020 were low, but not negligible, with 14.6M and 5.1M arrivals respectively in 2020 (Statistics Canada, 2021a). While essential workers were exempt from quarantine, all others were mandated a self-enforced 14-day quarantine upon re-entry. It was not until 22 February 2021 that a 3-day hotel quarantine was made mandatory for international arrivals by air (Government of Canada, 2021).”

p. 16, ll. 11-14: Ideally these nonstationary sampling probabilities by province (and really for all countries) would be built directly into the likelihood.

We agree with the reviewer that it would be interesting to incorporate sampling probabilities or proportions into either the likelihood of the internal node state or into the importation rates derived from them. We considered this, however decided against it for the following reasons. Pragmatically, the software used for ancestral character estimation is not flexible towards specifying priors on migration rates. We considered adjusting importation rates post hoc with the time-varying sampling probabilities by province and by country of origin, but doing so is subject to the same biases as our previous extrapolation of the upper estimate of the total number of introductions. We do not know the extent of travel bias or outbreak bias in the sequences, nor the case ascertainment rate, therefore we cannot simply assume that sequences were a random sample nor that importation rates scale linearly with available sequences. We will consider this suggestion for a future analysis using a software that can better accommodate extensions to the likelihood function.

Reviewer #2 (Recommendations for the authors):

1. Line 17-21 argue genomic epidemiology is useful to track and characterize VOI/VOC, but the sequences analyzed in this manuscript are largely pre-VOC so I think it makes sense to re-focus and expand this discussion to include commentary on the importance of evaluating the success of public health interventions, travel restrictions, etc.

We thank the reviewer for this recommendation and have changed the introduction to focus more on the precedence and importance of evaluating effectiveness of public health interventions.

“Characterization of viral importations over time can also clarify the effectiveness of public health interventions by associating inflection points in importations with drastic changes in policies such as international travel restrictions (Magalis et al., 2020). […] Evaluating travel restrictions’ effectiveness towards reducing burden could inform the stringency and timing of future public health interventions.”

The discussion also elaborated upon pragmatic policy changes that could be informed by these evaluations:

“The burden of cases in the second wave was primarily due to sublineages introduced in the summer of 2020, despite a low number of importations and case burden at that time. […] These advanced planning steps also serve to minimize the economic harm due to work absenteeism and supply chain disruption.”

2. I think outlining a hypothesis at the end of the intro is useful, but as it is currently written the hypothesis is very long and difficult to understand. Can the authors restate and clarify the hypothesis into a single sentence?

Thank you for the suggestion. We have distilled the scope into a single hypothesis:

“We tested the hypothesis that international travel restrictions enacted in March 2020 effectively reduced international importations of SARS-CoV-2 into Canada, yet ongoing introductions contributed to COVID-19 persistence into early 2021, exacerbated by highly transmissible B.1.1.7 and other VOC sublineages.”

3. As the authors discuss, their results and those of all genomic epidemiological studies are very sensitive to sequence quality, availability, distribution, and sampling paradigms. They reduce the impact of unequal global sequencing efforts by subsampling global sequences proportionally to their monthly case counts, excluding Canadian sequences. I understand the rationale here, but worry that inclusion of all available Canadian sequences (which are not proportional to case count trends as is clear in Figure 1D and 1E) skewed their results. Perhaps one way to assess the robustness of the results would be to repeat the main analyses using a subsampled (proportion to case counts) Canadian dataset and compare results? Relatedly, even case counts could be skewed depending on testing availability and accuracy -- could the authors include a brief commentary on the availability, accessibility, and capacity of testing centers in Canada during the first year of the pandemic?

We fully agree with the reviewer and have undertaken sensitivity analyses to address these concerns, as cited in text above. Additionally, in regards to the potential bias towards identifying USA importations across subsampling strategies, we added the following remark.

“Although the USA was highly represented in all of our subsamples as a result of its large contribution to COVID-19 cases in 2020 and high sequence availability during early months, our results suggest a greater effect than due to sampling alone. On average, the USA sequences represented 28.9% (28.7 – 29.2%) of total international sequences, yet accounted for 46.3% (44.0 – 48.7%) of all sublineages and 57.7% (55.6 – 59.8%) of singletons. Upon maximizing the number of Canadian sequences in the analysis, where global sequence representation was more normalized but less comprehensive, the USA sequences represented fewer of the international sequences (25.8%, 25.6 – 26.1%) and still accounted for 38.4% (37.0 – 39.8%) of sublineages and 46.4% (44.6 – 48.3%) of singletons.”

– The Alberta and Maritimes outbreaks in spring of 2020 appear very striking (large outbreaks with lots of person-to-person spread and very few descendants) in Figure. However, this interpretation is fraught when you consider Figure 1E and 1D because there are lots of available sequences during this time compared to the total case count vs fall of 2020 where there is actually a higher case burden and very few available sequences. I think many readers will get hung up on this so maybe the authors could address it briefly in the text?

We agree that the relative availability of sequences in the first and second waves is a limitation for several provinces, particularly Alberta. In the updated analysis, there was an increased representation of the second wave, improving this issue somewhat.

4. Lines 6-10. I think it is a big leap to estimate the total possible introductions based on sequences representative of 1.1% of total cases.

We agree with the reviewer that this upper estimate required further consideration of travel-related bias, outbreak-related bias, and case ascertainment rates. Since these values are not known across provinces, we have removed this upper estimate and replaced it with a discussion of the type of data needed to extrapolate the true number of introductions.

“Low sequence representation can lead to underestimates of total introductions if neither index case or descendants were sampled, underestimates of sublineage size if not all descendants were sampled, and similarly, overestimates of the proportion of singletons, which may have been from unsampled transmission chains. […] Additionally, prospective cohort studies or seroprevalence studies would ameliorate our estimate of the case ascertainment fraction.”

5. It might be very informative to take a more careful look at rates of viral spread (R0) over time in addition to patterns of viral importation. This could help contextualize the effect of public health interventions on domestic spread. The authors bring up quarantine, contact tracing, testing, and individual behavior change in the discussion -- do the authors think the rate of international imports is the best measure to assess these interventions? I think R0 would be just as informative!

We thank the reviewer for pointing out that showing differences in domestic viral spread over time would strengthen the conclusions regarding the effectiveness of NPIs. In regards to what is the best measure of the effectiveness of travel restrictions, we maintain that looking at importations over time is paramount. We experimented with modelling the effective reproduction number over time using national and provincial incidence, as well as sublineage descendants sampled over time, however there was considerable uncertainty in these estimates due to uncertainty in the serial interval and when occurrences were infrequent. We are developing more refined methods of deconvoluting the incidence data to improve these estimates for a follow-up publication. Instead, to address the reviewer’s request for a more thorough characterization of domestic spread, we have added a figure showing the distribution of sampled Canadian descendants of key sublineages (Figure 2B) and respective subtrees (Figure 2—figure supplement 1).

6. I am skeptical the results on reductions in sublineage size, increasing tMRC, and perhaps even sources of domestic transmission will persist if the Canadian sequences are subsampled proportionally to case counts. Can the authors please assess the robustness of these results against a subsampled dataset?

We have undertaken a Canada sequence subsampling sensitivity analysis and evaluated the effects on sublineage size, TMRCA, and relative sources of domestic transmission analysis. Please see subsampling sensitivity analysis.

7. I agree with and appreciate the authors' discussion of the importance data sharing, including a limited metadata set to aid in the rapid interpretation and use of sequence data for important public health decisions.

We thank the reviewer for sharing our opinion of how data sharing augments our collective knowledge of the epidemic and better positions us to respond effectively.

8. Although this is a stylistic preference, I believe a cleaned up version of Figure 7 with the Canada's key public health interventions and COVID-19 outbreaks along with a map of Canada's provinces (already a supplemental) and case counts during the study period would be very useful to include as an early, primary figure.

Thank you – we agreed with the reviewer’s recommendation to move the summary of interventions, stringency, cases, and sequences to be the first primary figure, as it gives an overview of the epidemic and sets the stage for the analysis. We also added the Canadian map in place of the legend in Figure 1A to provide geographic context for provincial names.

Reviewer #3 (Recommendations for the authors):

As mentioned in the Public Review, my main concern lies on the explicit naming of sources of importation following the application of an epidemiologically-based approach to downsample the publicly available genome sequences from GISAID. While it is mentioned that ten subsamples were generated, it is not clear to me that the ancestral state reconstruction was evaluated on these ten subsamples and how the results compare between them. The uncertainty of the ASR itself was incorporated for the fixed phylogeny through an ML approach (ape::ace()), but I am missing the analysis of the robustness of the results across different subsamples.

We thank the reviewer for noting that we could better show the uncertainty across bootstraps. All of the reported estimates in the Results are the mean across ten bootstraps. Including uncertainty in every figure was not feasible as there were already many features in the plot, which we did not want to detract attention from. In the subsampling sensitivity analysis undertaken, the uncertainty across bootstraps for each subsampling strategy was shown more explicitly.

Furthermore, even if these analyses where performed, it is unclear whether a different subsampling approach could lead to the identification of a different profile of source countries, and if so, which subsampling approach would be more accurate. The accuracy of the identified countries of origin (and therefore accuracy of the proposed downsampling approach) would be validated if a complementary data source showed similar patterns: an expected importation index (EII ~ Incidence * no. of travellers to Canada), self-reported travel histories of patients with confirmed COVID-19, etc. If these complementary data sources are unavailable, I would recommend reducing the emphasis on the identification of countries of origin (or constraining it to regions where you're confident that the monthly reported cases are reliable – the USA would be an obvious example of this), or adding a clear discussion on these limitations. I would argue that this goes beyond the already mentioned idea that downsampling does not allow to increase the numbers of sequences from some unsampled times/locations (page 16, lines 5-20), and that it can also reduce the probability of phylogenetically identifying true importations from countries with low case prevalence.

On other topics: I am unsure about the estimation of total number of introductions presented in page 8 (lines 6-10). Do you have any particular reason why you would expect that importations scale linearly with sampling/representation? Given the overall dynamics, wouldn't it be more probable that higher sampling intensities would add unsampled sequences to the larger sublineages in the country (proportional to the size of the lineages) as well as uncovering unsampled importations? The rationale behind these calculations requires some more detailed explanations.

We agree with the reviewer that the estimated total number of introductions was too simplistic and rife with assumptions. We have omitted this upper estimate and expanded upon why this is challenging to estimate in the discussion.

“Low sequence representation can lead to underestimates of total introductions if neither index case or descendants were sampled, underestimates of sublineage size if not all descendants were sampled, and similarly, overestimates of the proportion of singletons, which may have been from unsampled transmission chains. […] Additionally, prospective cohort studies or seroprevalence studies would ameliorate our estimate of the case ascertainment fraction.”

Regarding the domestic transmission, I was wondering if you had access to finer geographical resolution beyond province-level? Provinces are particularly large geographical areas, and providing some insights (rather than a full, high-resolution phylogeographic re-analysis of all the data) into the sampling within provinces could add to the manuscript. This is particularly the case for the analysis of international importations to Quebec and Ontario: where these sublineages first detected in large urban areas or not necessarily? Is the sampling within provinces representative of the cases reported in urban and rural areas? Are there any differential dynamics between larger urban areas and more remote locations?

Unfortunately, we did not have access to geography metadata for Canadian sequences at a finer scale than province. This information is not available in the public domain presumably for confidentiality reasons. Our group has collaborated with the BCCDC to model importations and transmission among the five health authorities in BC. We have also worked with Ontario public health colleagues to model relative contributions of international transmission using sequences available from over a dozen provincial health regions as part of an extensive sampling effort on 24 November 2020 for a point prevalence study. The manuscript for this work is still in preparation and we will take your suggestions into consideration therein, stratifying based on rural or urban geographic landscape. I would expect that generally urban centers act as major sources of viral transmission (based on density and travel proclivity) to more remote locations. However, we are not able to comment on this dynamic based on the data available in this study.

Finally, it is quite interesting to see that the peak of TMRCAs nicely coincide with the maximum stringency on March 21 (Figure 4A,B) – even more so given that the TMRCA of a sublineage is more likely representing something akin to the first transmission event within said sublineage rather than the importation event per se. It would be interesting if the TMRCAs could be contextualised to show their distribution after accounting for sampling lag (i.e., the time between symptom onset and sample collection). This could be estimated from sequence metadata or approximated from non-genomic epidemiological data, but it would give a better idea of the times during which the infections that lead to these sublineages started.

We thank the reviewer for distinguishing sublineage importation from TMRCAs. While we considered your suggestion to estimate and correct for sampling lag, we were hesitant to apply a fixed value or randomly draw from a distribution of lags in light of unknown variability described in the discussion below. Rather, we clarified the interpretation of sublineage TMRCAs.

“Sublineage TMRCAs typically lag importation by several days, but could precede them if multiple infections with subsequently sampled descendants occurred during travel. While accounting for sampling lag between symptoms and sample collection would get closer to the importation date, this assumes symptom onset occurred during travel, which might not always be the case. In contrast, if we assume that infected travellers who sparked domestic transmission chains were tested upon or near after arrival, then this lag would be zero to a couple days. If an infectious traveller was not sampled, but gave rise to an outbreak that was later sampled, then an additional lag must be considered regarding the length and number generations of unsampled viral predecessors since importation. A TMRCA could precede importation if the sublineage index case infected multiple people during travel. Where travel history is available, incorporating internal nodes or branch states for travel location would help to inform the timing as well as geographic origins of importations (Hong et al., 2021). In light of these additional sources of uncertainty, we abstained from applying an importation lag to TMRCAs, and emphasize their interpretation as the estimated time of the recent common ancestor of all sampled descendants.”

[Editors' note: further revisions were suggested prior to acceptance, as described below.]

The manuscript has been dramatically improved---the sensitivity analyses are immensely helpful---but there are some remaining issues that need to be addressed, as outlined below. The reviewing editor (also reviewer 1) is concerned that some of the policy suggestions do not follow the evidence presented.

We are highly grateful to the Senior Editor and Reviewing Editor for their comments and suggestions, which we have addressed below.

1. Abstract: "… suggesting travel restrictions and quarantine must be sustained to fully curtail COVID-19 burden." I am guessing most readers will likely interpret "to fully curtail COVID-19 burden" as to get the burden near zero. Is this what the authors intended? It affects the logic. The results here show that travel restrictions were usually associated with at least several importations per week. Mathematically, we know that the burden of COVID-19 does not really depend on the number of importations if R ~ 1 or more---all you need is at least one endemic lineage to maintain the prevalence/burden dictated by R. Maybe the authors are referring to conditions when R<1 and travel restrictions are even greater than what they have been historically, but this is unclear. It is also unclear precisely what the authors mean by "quarantine." Usually, it refers to the isolation of contacts of potential cases (not the cases themselves), and sometimes it refers more generally to lockdowns. Neither policy seems particularly worth calling out here. Why not say "sustained interventions to lower transmission (R<1)"? It seems unjustified for this study to single out quarantine (either def) as the necessary complementary intervention when so many other approaches have been shown more effective. I'm dwelling on the wording because these kinds of conclusions, especially in an abstract, can have an outsized impact on policy.

We thank the reviewer for identifying this statement as problematic. We have changed the wording as follows:

“Despite the drastic reduction in viral importations following travel restrictions, newly seeded sublineages in summer and fall 2020 contributed greatly to the persistence of COVID-19 cases in the second wave, highlighting the importance of sustained interventions to reduce transmission.”

2. Abstract: "…restrictions that reduce the probability of importations are most effective during periods of low domestic prevalence and low or waning immunity." Waning immunity means R is rising, which, all things equal reduces extinction probabilities and increases the success of invading lineages (e.g., Otto and Whitlock 1997). Without more precision about R and how effectiveness is being defined, it's not clear restrictions should be more impactful under these conditions. I do see the positive correlation between "periods of waning immunity" and possible ban-assisted extinction, but only because immunity is most often waning when there has recently been an epidemic, and the recently boosted immune protection is relatively high. I worry again about the implications of these claims, which are not exactly supported by the research here.

We agree that we have extrapolated our claims beyond what was supported by our findings, where we did not explicitly model changes in population-level immunity. This sentence has been modified accordingly:

“Although viral importations are nearly inevitable when global prevalence is high, with fewer importations there are fewer opportunities for novel variants to spark outbreaks in new communities or outcompete previously circulating lineages.”

3. Relatedly, the main text describes evaluating the relative contribution of (reducing) importations to (reducing) *burden* (e.g., ll. 190-192, ll. 605-607, 691-701), but this analysis is not really performed here. Such analysis would require careful estimation of R and prevalence and a good enough transmission model to evaluate extinction probabilities and outbreak sizes (in the case of R<1) over time, and it would have to account for clonal interference. Put differently, the fact that imported lineages caused x% of cases does not mean that the case count would've been (100-x)% without those importations. The fact that imported lineages displaced resident lineages does not mean the resident lineages could not have caused comparable infections had the imported lineages been stopped at the border, assuming the imported lineages had no fitness differences. Imperfect travel bans that allow some lineages to invade can at best slightly (at measured rates) change the timing of an epidemic (ll. 704-707), but usually not its size or severity. The authors write, "Although blind travel bans may not be beneficial from a socioeconomic perspective, ultimately governments need to protect citizens and dynamic travel bans are one of the few tools available." This is surprising considering the number of broader tools at our disposal (e.g., rapid testing, vaccination, ventilation/filtration) that can slow importations AND endemic spread, as well as the observed ineffectiveness of travel bans at stopping spread outside of islands and China. The authors IMO have not demonstrated how much bans could delay importation and allow time to "plan", or if investments to reduce transmission generally are more effective. I suggest these claims be carefully reconsidered.

I apologize for the wordiness. I would be happy to bolster these points (maybe with math) if they are unclear. I am concerned that too much is being extrapolated about the possible and actual epidemiological impacts of travel bans from the measured rates of importation.

We appreciate the reviewer’s perspective on the semantics and interpretation of our findings. We have not asserted that the case count correlates directly with importations, but that each importation increases the chance of seeding a new community and the chance of introducing a more transmissible variant. Wording has been adjusted as follows:

ll. 190-192: “Evaluating travel restrictions’ effectiveness towards reducing importations could inform the stringency and timing of future public health interventions.”

ll. 605-607: “A more rapid and stringent public health response would have reduced the initial burden by preventing early sublineages from establishing widespread transmission chains that persisted throughout the first and second waves.”

ll. 691-705: “Case burden in the second wave was primarily due to sublineages introduced in the summer of 2020 amid low prevalence, despite a low number of importations. This suggests that ongoing seeding events impacted viral persistence prior to the arrival of VOCs with increased transmissibility. Considering the global and temporal diversity of the SARS-CoV-2 pandemic, it is difficult to predict what variant may be introduced and subsequently disperse through a local population. While broad and longstanding restrictions against non-essential international travel is not necessarily an advisable policy in light of economic impacts, swift and stringent travel bans towards locations harbouring a high frequency of a VOC not yet identified domestically should be seriously considered to reduce the probability of multiple, simultaneous outbreaks being seeded and potentially overwhelming the healthcare system.”

ll. 704-707: “Dynamic travel bans can be detrimental from a socioeconomic perspective; therefore their implementation must be weighed relative to other non-pharmaceutical interventions to reduce transmission.”

You may wish to consider some suggestions to improve readability from another reviewer:

4. In line 424, I'm not entirely sure what the authors mean by "the first wave of singletons had a higher amplitude than the second wave". By 'amplitude', do they refer to the range of dates when singletons were identified between both waves? Is this relative to a feature of the importations of sublineages (given that this sentence talks about this comparison between sublineage and singleton dynamics)? Some clarification would be helpful.

Thank you for seeking clarification on the wording here. By amplitude, we meant the wave height of the weekly singletons introduction rate, or the total average number of singletons introduced per week. Wording changed as follows:

“Importation dynamics for singletons mostly mirrored sublineage trends (Figure 4—figure supplement 1; Figure 3—figure supplement 1), although the USA contributed relatively more towards singletons than sublineages in both waves, and the maximum singleton importation rate was higher in the first wave than the second wave.”

5. The Discussion addresses a lot of the reviewer comments thoroughly and systematically, which was very helpful, but could be streamlined for the readers. I think that, for example, the two paragraphs that go from line 711 to line 741 could be reduced to a single paragraph discussing 1. sampling bias in phylogeography, 2. the specific effects of overrepresentation of some locations and how they addressed it with subsampling and the sensitivity analyses, 3. effects of upsampled lineages and under-sampled locations and how this can theoretically be overcome, 4. the broader challenges of scaling up Bayesian phylogenetics (this last point could be a separate paragraph, but unrelated to sampling bias).

We are grateful for this feedback towards improving the flow of the Discussion and have modified the text as follows:

“Ancestral trait reconstruction is sensitive to sampling bias as ancestral node states are more likely to reflect states predominant among tips (Baele et al., 2016; Hill et al., 2021; Lemey et al., 2009). […] Stratified analyses of well-supported (sub)lineages in future analyses will facilitate direct comparisons of Bayesian and maximum likelihood phylogeography methodologies towards estimating importation rates.”

Author details

1. Angela McLaughlin

   1. British Columbia Centre for Excellence in HIV/AIDS, Vancouver, Canada
   2. Bioinformatics, University of British Columbia, Vancouver, Canada

   Contribution

   Conceptualization, Data curation, Software, Formal analysis, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing – original draft, Project administration, Writing – review and editing, AM conceived of the project, conducted most of the data cleaning, analysis, and visualization, and wrote the original draft of the manuscript

   For correspondence

   amclaughlin@bccfe.ca

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0001-5606-9080

2. Vincent Montoya

   British Columbia Centre for Excellence in HIV/AIDS, Vancouver, Canada

   Contribution

   Data curation, Formal analysis, Methodology, Writing – review and editing, VM helped to develop methodologies for data cleaning, alignment, and tree building

   Competing interests

   No competing interests declared

3. Rachel L Miller

   1. British Columbia Centre for Excellence in HIV/AIDS, Vancouver, Canada
   2. Bioinformatics, University of British Columbia, Vancouver, Canada

   Contribution

   Visualization, Writing – review and editing, RLM contributed towards improving the graphics

   Competing interests

   No competing interests declared

4. Gideon J Mordecai

   Department of Medicine, University of British Columbia, Vancouver, Canada

   Contribution

   Visualization, Methodology, Writing – review and editing, GJM provided feedback on the phylogeographic methods and graphics, as well as provided citations

   Competing interests

   No competing interests declared

5. Canadian COVID-19 Genomics Network (CanCOGen) Consortium

   Contribution

   Funding acquisition, Writing – review and editing, Initial feedback on earlier versions of the analysis

   Competing interests

   No competing interests declared

   Additional information

   See supplementary materials for list of consortium members and affiliations.

6. Michael Worobey

   Department of Ecology and Evolution, University of Arizona, Tucson, United States

   Contribution

   Conceptualization, Supervision, Funding acquisition, Visualization, Methodology, Writing – review and editing, MW contributed suggestions on estimating uncertainty, visualizations, and interpretations

   Competing interests

   No competing interests declared

7. Art FY Poon

   Department of Pathology and Laboratory Medicine, Western University, London, Canada

   Contribution

   Data curation, Supervision, Visualization, Methodology, Writing – review and editing

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0003-3779-154X

8. Jeffrey B Joy

   1. British Columbia Centre for Excellence in HIV/AIDS, Vancouver, Canada
   2. Bioinformatics, University of British Columbia, Vancouver, Canada
   3. Department of Medicine, University of British Columbia, Vancouver, Canada

   Contribution

   Conceptualization, Resources, Data curation, Software, Supervision, Funding acquisition, Validation, Visualization, Methodology, Writing – original draft, Project administration, Writing – review and editing, JBJ contributed to refining the research questions

   For correspondence

   jjoy@bccfe.ca

   Competing interests

   No competing interests declared

   "This ORCID iD identifies the author of this article:" 0000-0002-7013-1482

• 1,855

  Page views

• 301

  Downloads

• 4

  Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.