Regional importation and asymmetric within-country spread of SARS-CoV-2 variants of concern in the Netherlands

  1. Alvin X Han  Is a corresponding author
  2. Eva Kozanli
  3. Jelle Koopsen
  4. Harry Vennema
  5. RIVM COVID-19 molecular epidemiology group
  6. Karim Hajji
  7. Annelies Kroneman
  8. Ivo van Walle
  9. Don Klinkenberg
  10. Jacco Wallinga
  11. Colin A Russell
  12. Dirk Eggink
  13. Chantal Reusken  Is a corresponding author
  1. Department of Medical Microbiology and Infection Prevention, Amsterdam University Medical Centers, Netherlands
  2. Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Netherlands

Abstract

Background:

Variants of concern (VOCs) of SARS-CoV-2 have caused resurging waves of infections worldwide. In the Netherlands, the Alpha, Beta, Gamma, and Delta VOCs circulated widely between September 2020 and August 2021. We sought to elucidate how various control measures, including targeted flight restrictions, had impacted the introduction and spread of these VOCs in the Netherlands.

Methods:

We performed phylogenetic analyses on 39,844 SARS-CoV-2 genomes collected under the Dutch national surveillance program.

Results:

We found that all four VOCs were introduced before targeted flight restrictions were imposed on countries where the VOCs first emerged. Importantly, foreign introductions, predominantly from other European countries, continued during these restrictions. After their respective introductions into the Netherlands, the Alpha and Delta VOCs largely circulated within more populous regions of the country with international connections before asymmetric bidirectional transmissions occurred with the rest of the country and the VOC became the dominant circulating lineage.

Conclusions:

Our findings show that flight restrictions had limited effectiveness in deterring VOC introductions due to the strength of regional land travel importation risks. As countries consider scaling down SARS-CoV-2 surveillance efforts in the post-crisis phase of the pandemic, our results highlight that robust surveillance in regions of early spread is important for providing timely information for variant detection and outbreak control.

Funding:

None.

Editor's evaluation

Han et al., analyze sequences from randomly sampled COVID-19 cases in the Netherlands to understand the impact of flight restrictions on the importation of SARS-CoV-2 variants. In line with prior observations and common wisdom, they find that targeted flight restrictions were not effective at preventing introductions of new lineages and that their early spread in the Netherlands was sustained by urban centers. These useful findings, based on unusually strong sequence collection techniques, can inform surveillance policy and improve basic understanding of the spread of SARS-CoV-2 variants.

https://doi.org/10.7554/eLife.78770.sa0

Introduction

Coronavirus disease 2019 (COVID-19) has resulted in excess morbidity and mortality across the world. In response, governments have implemented travel restrictions and nonpharmaceutical interventions in order to limit introductions and reduce transmission of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2; Brauner et al., 2021; Chinazzi et al., 2020; Flaxman et al., 2020). However, high levels of global infections have led to the evolution and emergence of variants of concern (VOCs) that are more transmissible, some of which encode putative mutations that evade immunity acquired from previous infection or vaccination (Harvey et al., 2021). These VOCs have led to resurging SARS-CoV-2 outbreaks, hampering efforts to contain and control the pandemic worldwide. Of note, four such VOCs arose into global prominence in late 2020, including Alpha (Nextclade 20I; PANGO lineage B.1.1.7), Beta (20 H; B.1.351), Gamma (20 J; P.1), and Delta (21 J; B.1.617.2), causing substantial levels of transmission worldwide, with Alpha and Delta being the most common VOCs globally in 2021 (Campbell et al., 2021).

The Alpha VOC was first reported in the United Kingdom (U.K.) during the fall of 2020 and found to be 43–90% more transmissible (Davies et al., 2021a; Volz et al., 2021) with greater mortality risks (Davies et al., 2021b; Grint et al., 2021) than previously existing variants. Of the 17 amino acids mutations found in Alpha, N501Y in the receptor-binding domain (RBD) of the spike protein was predicted to increase binding to the human angiotensin-converting enzyme 2 receptors (Starr et al., 2020). This is also a common mutation found in Beta (Tegally et al., 2021) and Gamma (Faria et al., 2021). On the other hand, the Delta VOC, first identified in India in October 2020 (Cherian et al., 2021), encodes P681R mutation in the furin cleavage site in spike protein and R203M mutation in the nucleocapsid protein that improves infectivity (Syed et al., 2021). Delta has also been linked to increased disease severity, as well as greater and longer viral shedding (Ong et al., 2021). In the U.K., where the VOC was first detected in April 2021, epidemiological modelling estimated the VOC to be 40–80% more transmissible than Alpha.

The four aforementioned VOCs were also introduced into the Netherlands, with the Alpha and Delta VOCs subsequently dominating infections in the country in 2021. In a bid to deter introductions and slow down the spread of VOCs, the Dutch government implemented targeted flight restrictions on countries where these VOCs had first emerged. Various non-pharmaceutical interventions were also implemented as the country experienced multiple waves of infections between 2020 and 2021. Since the end of 2020, the Dutch National Institute for Public Health and Environment scaled up its sequencing efforts under a random national surveillance program. This detailed surveillance program allows the monitoring of the introduction and spread of novel variants or specific mutations. Here, SARS-CoV-2 positive samples were randomly collected across the country and 39,844 high-quality SARS-CoV-2 whole genomes were sequenced between 22 September 2020 and 31 August 2021 (48 calendar weeks). We analyzed these sequences alongside epidemiological data to identify importation risks of novel variants into the Netherlands and characterize their subsequent patterns of spread within the country.

Methods

Whole genome sequencing

A total of 39,844 nasopharyngeal samples were randomly collected across all 25 GGD municipal health services across all municipalities in the Netherlands and were sequenced for whole SARS-CoV-2 genomes. These samples were collected through nationwide community testing programs including at test buses and health facilities (https://www.rivm.nl/coronavirus-covid-19/onderzoek/kiemsurveillance). Test-positive samples were randomly subsampled each week in each municipality to minimize case ascertainment bias. Only specimens with cycle threshold (Ct) values <30 were selected for whole genome sequencing.

As testing and samples were analyzed by 30–35 different laboratories across the country, different nucleic acid extraction methods were used (Herrebrugh et al., 2021). For samples analyzed by the laboratory of the Dutch National Institute for Public Health and Environment (38,260 samples; 96% of all samples analyzed), total nucleic acid was extracted using MagNApure 96 (MP96) with total nucleic acid kit small volume (Roche). RT-qPCR was performed on 5 μl total nucleic acid using TaqMan Fast Virus 1-Step Master Mix (Thermo Fisher) on Roche LC480 II thermal cycler with SARS-like beta coronavirus (Sarbeco) specific E-gene primers and probe and EAV as described previously (Corman et al., 2020; Scheltinga et al., 2005).

Amplicon-based SARS-CoV-2 sequencing for was performed using the Nanopore protocol ‘PCR tiling of COVID-19 virus (Version: PTC_9096_v109_revE_06FEB2020)’ which is based on the ARTIC v3 amplicon sequencing protocol (Artic Network, 2020). Several modifications were made to the protocol as primer concentrations were increased from 0.125 to 1 pmol for the following amplicon primer pairs: in pool A amplicons 5, 9, 17, 23, 55, 67, 71, 91, 97 and in pool B amplicons 24, 26, 54, 64, 66, 70, 74, 86, 92, 98. Both libraries were generated using native barcode kits from Nanopore SQK-LSK109 (EXP-NBD104, EXP-NBD114 and EXP-NBD196) and sequencing was performed on a R9.4.1 flow cell multiplexing 2 up to 96 samples per sequence run. Consensus sequences (>50 x depth-of-coverage) are generated using an in-house bioinformatic pipeline developed by the Dutch National Institute for Public Health and Environment (https://github.com/RIVM-bioinformatics/SARS2seq/; Zwagemaker, 2022).

Epidemiological data

All epidemiological data including the breakdown of positive cases by age group and weekly number of laboratory-confirmed cases in each Municipal and Regional Health Service region are provided by the Dutch National Institute for Public Health and Environment (https://www.rivm.nl/en/node/163991).

Phylogenetic analyses

We performed ancestral reconstruction analyses for each VOC lineage to identify likely overseas introduction into the Netherlands at the continental level, differentiating the Netherlands from the rest of Europe. As proportions of cases for each VOC lineage are unknown for most countries, we subsampled global sequences downloaded from GISAID (https://www.gisaid.org; dataset up to 6 October 2021) by the proportion of COVID-19 cases reported per week for each country using data from the Johns Hopkins University, Center for Systems Science and Engineering (CSSE) (http://github.com/CSSEGISandData/COVID-19). We aimed to sample 100 sequences each week. To ensure that countries that are underreporting cases (Gill et al., 2022) were included in our analyses, at least one representative sequence was included for each country with reported cases that week. The actual number of sequences sampled each week for each variant may differ because of sequence availability and the requirement to have at least one sequence from each country each week. For sequences from the Netherlands, we also subsampled Dutch sequences by randomly subsampling the number of sequences in each GGD region each week to the corresponding relative number of reported cases in the same week for that GGD region. As there may be certain weeks when our aforementioned strategy sampled less non-Netherlands sequences than Netherlands sequences, we would resample a larger number of non-Netherlands sequences accordingly such that there are at least twice as many non-Netherlands sequences as Netherlands sequences each week. The subsampling procedure yielded 6,365 (2,369), 1,531 (90), 1,274 (102) and 6,929 (1,035) Alpha, Beta, Gamma, and Delta global (Dutch) sequences, respectively.

All sequences were aligned to hCoV-19/Wuhan/WIV04/2019 (WIV-04; EPI_ISL_402124) using MAFFT v7.427 (Katoh and Standley, 2013). Likely problematic sites (https://github.com/W-L/ProblematicSites_SARS-CoV2) along with untranslated regions in the 5’ and 3’ ends were masked. All maximum-likelihood (ML) phylogenetic trees were reconstructed using IQ-TREE (Nguyen et al., 2015) under the Hasegawa–Kishono–Yano nucleotide substitution model with a gamma-distributed rate variation among sites (HKY +G). We regressed the root-to-tip genetic distances against sampling dates using treetime v0.8.1 (Sagulenko et al., 2018) to assess the level of temporal signal, ensuring that none of the representative sequences were deemed molecular clock outliers.

Using these sequences, we then reconstructed approximate ML phylogenies using FastTree v2.1.11 (Price et al., 2010). All phylogenetically neighboring overseas sequences placed within two nodes away from any Dutch sequence were retained. This further reduced the number of sequences to a representative set of 3671, 496 and 575 and 2180 sequences for Alpha, Beta, Gamma and Delta VOCs respectively. We then reconstructed an ML phylogenetic tree again under the HKY +G nucleotide substitution model using IQ-TREE after removing any molecular clock outlying sequences identified by treetime. Here, however, we included the WIV-04 reference genome in the phylogeny reconstruction which was used as an outgroup for tree rooting. We then time-scaled these ML phylogenies using treetime, which were then used as fixed tree topologies in BEAST v.1.10.4 (Suchard et al., 2018) to perform Bayesian discrete phylogeographical analyses at the continental level. Here, we performed 100 million MCMC generations, sampling every 1000 steps.

We randomly downsampled sets of Alpha (n=1389) and Delta (n=1342) variant sequences collected in the Netherlands. To understand within-country source-sink dynamics during early introductions and proliferation patterns during later periods, we used BEAST v.1.10.4 to perform continuous phylogeographical analyses on these sequence data, using a relaxed random walk diffusion model and a Cauchy distribution model among branch heterogeneity in diffusion velocity (Lemey et al., 2010). We used HKY + G nucleotide substitution model and a Skygrid coalescent model (Gill et al., 2013; each grid point denoting one week) with a relaxed lognormal molecular clock. We inferred geographical coordinate input using the first four digits of postcodes (i.e. neighborhood level) associated with the sampled sequences. For sequences with identical postcodes, we randomly selected geographical coordinates corresponding to the postcode area using shapefiles provided by https://www.gadm.org. We performed two independent chains, each with 300 million MCMC generations for each variant analysis, sampling every 50,000 steps. The first 100 million steps were discarded as burn-in. Assessment of convergence (effective sample size >200) for each chain was performed using Tracer v1.71 (Rambaut et al., 2018). The posterior distribution of trees from both chains were combined and resampled at every 100,000 steps to generate the maximum clade credibility tree. Customized scripts from the SERAPHIM package (Dellicour et al., 2016) were used to extract the inferred phylogeographic reconstruction.

All tree visualizations were performed using baltic (https://github.com/evogytis/baltic; Dudas, 2021).

Aggregated mobility data

We used publicly available mobility data from Google COVID-19 community mobility reports (https://www.google.com/covid19/mobility/) which contain daily anonymized location histories as a measure of people’s movements. Google mobility data consisted of six categories that were measured relative to a baseline value. This baseline is the median mobility value between pre-pandemic weeks of 3 January and 6 February 2020. Categories include residence, parks, retail and recreation, groceries and pharmacies, working place and transit. Data for different regions of the Netherlands were available. We calculated aggregated nationwide mean mobility by averaging values across all regions for all categories except for residence and parks where the former has a reversed effect on relative mobility while the latter is affected by climate.

Results

SARS-CoV-2 infections and genotypes circulating in the Netherlands from September 2020 to August 2021

There were 1,792,759 laboratory-confirmed SARS-CoV-2 cases in the Netherlands during the study period between 22 September 2020 and 31 August 2021 (week 39/2020 to week 34/2021; Figure 1A). Similar to the first wave of the pandemic in the Netherlands in Spring 2020, most reported cases were attributed to the more densely populated regions of the country including North and South Holland, as well as North Brabant (Figure 1B) where the first local clusters of SARS-CoV-2 were also detected in March 2020 (Oude Munnink et al., 2020). 39,844 SARS-CoV-2 positive nasopharyngeal samples were randomly selected from 25 Municipal and Regional Health Service (GGD) regions across the Netherlands during this study period and sequenced to obtain whole virus genomes as part of the national SARS-CoV-2 genomic surveillance program.

SARS-CoV-2 infections and VOCs circulating in the Netherlands from September 2020 to August 2021.

(A) Weekly number of laboratory-confirmed SARS-CoV-2 cases (1st panel from the top) and sequenced genomes (2nd panel). Lineage proportions of sequences colored by NextClade genotype designations (3rd panel). Breakdown of positive cases by age group from data provided by the Dutch National Institute for Public Health and Environment (4th panel). Aggregated weekly average percentage change in mobility to the baseline in the Netherlands from Google’s COVID-19 community mobility reports. Baseline mobility is the median value from a 5-week period between 3 January 2020 and 6 February 2020, prior to the COVID-19 pandemic in Europe (5th panel). (B) Mean number of laboratory-confirmed cases per 100,000 inhabitants (data from the Dutch National Institute for Public Health and Environment; left panel); total number of sequenced genomes in different Municipal and Regional Health Service (GGD) regions over the entire study period (middle panel); Percentage of Dutch population residing in each GGD region (right panel).

Using NextClade lineage assignment (Aksamentov et al., 2021), the viruses sampled at the start of the study period were largely genotyped as clade 20 A and its daughter lineages, 20B and 20E (EU1) (Figure 1B–C). 20 A was the lineage that seeded the pandemic in Europe in March 2020. On the other hand, 20E (EU1) was first detected in Spain on June 2020 and spread widely across Europe due to the resumption of regional travel over summer 2020 (Hodcroft et al., 2021). Owing to rising case numbers, non-pharmaceutical interventions closing restaurants and nightlife establishments were implemented on 14 October 2020. Cases dipped momentarily while both 20 A and 20E (EU1) remained co-circulating into December 2020.

The first Alpha sample was collected on 5 December 2020 in the national surveillance program prior to the full lockdown that closed all public venues, workplaces and schools on 15 December 2020. A curfew was also imposed later on 23 January 2021. A sharp drop in cases was observed after the implementation of the full lockdown. Alpha then displaced 20 A and 20E (EU1) over time to become the dominant circulating virus lineage by 16 February 2021 (week 6) for the rest of the lockdown period. Other VOCs such as Beta (N=422 sequences; first sequence was collected on 22 December 2020) and Gamma (N=350 sequences; first sequence was collected on 27 January 2021) were also detected by random surveillance in The Netherlands around the turn of the year but did not circulate to the levels of Alpha.

Alpha caused a rebound in cases around mid-March 2021, after which case numbers stabilized and eventually began to decline at the end of April 2021. The Dutch government began taking steps to relax restrictions around the same time, starting with the end of curfew and resumption of higher education during the week of 27 April 2021. The first Delta sample was collected in the previous week on 15 April 2021 and continued to accumulate in frequency. By week 25 (29 June 2021), the Delta VOC accounted for 24% of all weekly genomes sequenced.

Most restrictions were lifted in the same week by 26 June 2021. As SARS-CoV-2 prevalence declined over time, average nationwide mobility also increased steadily which eventually came close to pre-pandemic levels in June 2021 (mean percentage change relative to pre-pandemic baseline = –5.0% (s.d.=11.0%); Figure 1A). SARS-CoV-2 prevalence was at its lowest then with only 8,690 reported positive cases that week. Within only 1 week after reopening, however, weekly cases soared above 50,000 on weeks 26 and 27 (6–20 July 2021). With most infections attributed to Delta, the novel VOC replaced the Alpha VOC as the dominant lineage within the next 3 weeks as over 90% of randomly surveilled genomes were typed as Delta VOCs by mid-July.

Stratifying the number of weekly reported positive cases by patient age group, the relative proportions in case positive rates remained fairly consistent throughout the study period except for weeks 26 and 27 where the rapid increase in cases was largely attributed to individuals aged between 15 and 30 years (Figure 1A). One of the reasons behind widespread transmission among young adults then was super-spreading linked to nightlife venues (Koopsen et al., 2022). In response, the government shut nightlife establishments down again on 10 July 2021 (week 27). Case numbers fell promptly after but remained at over 30,000 new cases per week for the rest of the study period. The Delta VOC had in principle completely displaced Alpha by then with over 99% of randomly surveilled genomes sampled from August 2021 onwards.

Overseas introduction of variants of concern

To deter the introductions of novel VOCs into the Netherlands, travel restrictions were imposed on countries where the VOCs first emerged, including the United Kingdom between December 2020 and March 2021 due to the emergence of Alpha; South Africa and Brazil between January and June 2021 due to Beta and Gamma respectively; and India from April to June 2021 due to Delta. These travel restrictions include a ban on all incoming passenger flights except for those carrying cargo and medical personnel, on top of an entry ban for all non-European Union residents (Government of the Netherlands, 2021). On the other hand, travel within parts of the Schengen Area in Europe, which includes the Netherlands, remained largely possible during this period. To identify likely where and when VOCs were actually introduced into the country, we subsampled a representative set of Dutch and overseas sequences collected over the same time period. We then reconstructed time-scaled, maximum likelihood (ML) phylogenies and used these fixed trees to perform discrete trait analyses using a Bayesian approach to infer likely overseas introductions at the continental level. This was done by identifying subtrees subtending Dutch sequences with ancestral states that were attributed to an overseas origin (Figure 2). All four VOCs were already introduced into the Netherlands prior to the targeted flight restrictions that were imposed on countries where these VOCs first emerged.

Figure 2 with 3 supplements see all
Likely overseas introduction of VOC lineages into the Netherlands at the continental level.

For each VOC lineage, a time-scaled maximum likelihood phylogeny using the Dutch and their nearest overseas neighboring sequences was inferred. Discrete trait analyses were performed to infer the likely continental region of ancestral states. Subtrees or singletons with ancestral nodes attributed to an overseas origin but subtend only Dutch sequences are drawn. Shaded plot area denotes the timespan when a targeted flight restriction was imposed on the country where the VOC lineage first emerged (i.e. (A) Alpha, United Kingdom.; (B) Beta, South Africa; (C) Gamma, Brazil; (D) Delta, India).

Importantly, besides countries where travel restrictions were in place, we estimated multiple likely introduction events from other foreign countries into the Netherlands for all four VOCs (Alpha, n=100; Beta, n=7; Gamma, n=12; Delta, n=213). Given disparities in global sequencing efforts (Brito et al., 2021), the random surveillance strategy used in local sample collection, and low genetic diversity among SARS-CoV-2 genomes used to reconstruct ancestral states, we are unable to fully and reliably quantify the number of introductions attributed to different geographical regions. However, many of the estimated regions for these ancestral states were in Europe (Alpha, 71% of all estimated overseas introduction events; Beta, 29%; Gamma, 71%; Delta, 79%; Figure 2). Furthermore, these European introductions continue to occur during the targeted travel ban period. Inspecting the nearest phylogenetic ancestral taxon to the aforementioned subtrees, we found that many of these nearest overseas neighboring tips were detected in Belgium, Germany, France and Denmark where borders between the Netherlands remained open as well as other countries (e.g. Spain, U.K., Poland, U.S.) where no targeted travel restrictions were set in place (Figure 2—figure supplement 1). There was also no isolated period in time in which these VOCs were introduced into the Netherlands - introductions likely occurred repeatedly during the period when these VOCs were also proliferating within the country.

Within-country transmission dynamics of the Alpha and Delta VOCs

To further elucidate the transmission dynamics of the Alpha and Delta VOCs within the Netherlands, we performed continuous phylogeographic analyses using separate downsampled sets of Alpha and Delta sequence data (Figure 3). For the first four weeks since the initial detection of both VOCs within the country, introductions and phylogenetic branch movements were mostly concentrated in the more populous regions of the country, including North and South Holland, Utrecht and North-Brabant, forming a core of early dominant locations. During this period, dispersal events to regions outside of these GGD regions occurred as well but are relatively less frequent. However, as local infections were seeded in these areas, bidirectional exchanges in phylogenetic branches between different regions emerged throughout the country. These bidirectional exchanges continued to increase as prevalence of the VOC grew over time, even amidst a strict lockdown in the case of Alpha (Figure 4). In particular for the Delta VOC, we observed a rapid spike in inter-regional spread upon the week of 22–28 June 2021, with >400% estimated increase in total phylogenetic branch movements by 6 July 2021 (week 26). This significant rise in inter-regional exportation events likely contributed to the soaring case numbers observed between weeks 25 and 27 (22 June – 13 July 2021).

Source-sink dynamics the Alpha and Delta variants of concern (VOC) within the Netherlands during the first four weeks after their respective detection.

(A) Alpha (between 2 and 29 December 2020); (B) Delta (between 20 April and 18 May 2021).

Estimated number of phylogenetic branch movements and growth rate estimations of Alpha and Delta VOCs in the Netherlands.

Solid line shows the number of branch movements from early dominant source regions including North and South Holland, Utrecht and Brabant. Dashed line shows number of branch movements from areas outside of these early dominant source locations.

Discussion

Even if international travel restrictions are in place, the Netherlands is still highly vulnerable to importation risks of novel SARS-CoV-2 variants from its regional neighbors due to border policies within the European Union. As such, this regional vulnerability is not unique to the Netherlands and has been reported in other European countries as well (Lemey et al., 2021; Michaelsen et al., 2021; Osnes et al., 2021). Importantly, regional introductions of novel lineages often drive new waves of infections in Europe (Lemey et al., 2021). Prior to September 2020, the dominant variant lineages (i.e. 20 A and 20E (EU1)) that circulated the Netherlands were already seeded by imports from its European neighbours (Hodcroft et al., 2021; Nadeau et al., 2021; Oude Munnink et al., 2020). In fact, the initial introduction of SARS-CoV-2 in the Netherlands in February 2020 were attributed to travelers who visited Northern Italy where the earliest sustained European SARS-CoV-2 transmission network was seeded (Oude Munnink et al., 2020; Worobey et al., 2020). Here, we showed that all four VOC lineages detected in the country up to August 2021 also originated mainly from its European neighbors. Importantly, regional importation risks persisted throughout the period these VOCs circulated the country and overlapped with periods where targeted flight restrictions were imposed on countries where these VOCs first emerged. The emergence of the Omicron VOC in southern Africa in November 2021 (Viana et al., 2022) again led to reactionary targeted flight restrictions by several countries in the Global North, including the Netherlands which was still amidst a surging Delta infection wave. However, the Omicron VOC was already detected in samples collected one week before the imposed flight restriction and did not prevent it to rapidly become the dominant VOC circulating in the Netherlands by the end of 2021 (https://www.rivm.nl/coronavirus-covid-19/virus/varianten/omikronvariant). Previous studies showed that travel restrictions are only useful if restrictions barred arrivals from most countries provided that local incidence is low in the first place (Russell et al., 2021).

Due to disparities in global genomic surveillance efforts and the lack of travel history information among the sampled Dutch individuals, we could not make more precise and accurate phylogeographic inferences on overseas introduction of VOCs into the Netherlands (Lemey et al., 2020). Furthermore, there were countries that are underrepresented or missing in our subsampling of non-Netherlands sequences. While we ensured that at least one sequence from each country with confirmed cases and genomic data was included in our analyses (see Methods), there may be overseas introductions that we could not account for. We repeated our overseas phylogeographic analyses with an independent bootstrap subsample of sequence data (Figure 2—figure supplements 23). There are differences in the estimated distributions of countries with nearest overseas neighboring sequences to Dutch subtrees (Figure 2—figure supplements 1 and 3). However, our conclusions that there were multiple introductions from other European countries before and after travel restrictions were imposed on countries where VOCs first emerged remain evident in the bootstrap analysis (Figure 2—figure supplement 2).

We also found that early introductions of VOCs, specifically the Alpha and Delta VOCs, are more likely found in populous regions of the Netherlands, including Utrecht, North, and South Holland where larger cities are located that are also international and regional travel hubs. These areas constitute a core cluster of dominant source locations that also exported infections to the rest of the country during first few weeks after the VOC’s introduction into the Netherlands. As the number of infections in areas outside of these dominant source locations increase over time, bidirectional exchanges would also become more frequent. This type of asymmetric spatial spread dynamics had been previously shown in the U.K. as well and was found to enhance the intrinsic transmissibility of Alpha (Kraemer et al., 2021). Additionally, enhanced mobility has also been previously linked to the resurgence of outbreaks across Europe (Hodcroft et al., 2021; Lemey et al., 2021). Recent work also showed that increased mobility and population mixing drove the rapid dissemination of Delta in the U.K. (McCrone et al., 2021). While our analyses do not provide a causal relationship between the relaxation of non-pharmaceutical interventions and frequency of export events, the asymmetric exportation frequencies from dominant source locations, increased human mobility in weeks 25–27 across the country as well as the intrinsic higher transmissibility of Delta relative to Alpha likely all contributed to the widespread spike in cases in the Netherlands.

Novel and fitter variants of SARS-CoV-2 will likely continue to emerge in the future. Our results, along with others, show that unless well-coordinated actions are taken across Europe to mitigate importation risks (Ruktanonchai et al., 2020), targeted travel restrictions implemented by individual European countries will neither prevent nor slow down the introduction of novel variants. Our work also shows that early within-country spread of VOCs may be taken into future consideration in future genomic surveillance strategies, especially as countries are gradually considering scaling down SARS-CoV-2 surveillance efforts. Both the Alpha and Delta VOCs were first detected in the early dominant source locations, usually those that are more populous with greater international connections, and circulated mostly within these areas during the initial period after introduction. As such, a robust level of surveillance efforts should still be maintained in these dominant source locations to provide timely actionable information on novel variant detection as well as infection control. These surveillance efforts should encompass a minimal level of clinical diagnostic testing capacity be maintained to ensure clinical genomic surveillance remains sensitive enough for early detection of novel variants (Han et al., 2022). Wastewater surveillance could also be included to facilitate early variant detection and identify cryptic transmissions amid falling testing rates (Karthikeyan et al., 2022).

Data availability

All sequencing data have been deposited in the GISAID database (https://www.gisaid.org) under the accession codes listed in . All codes for our analyses are available at https://github.com/AMC-LAEB/nl_sars-cov-2_genomic_epi_2022; Han, 2022.

Data availability

All sequencing data have been deposited in the GISAID database (https://www.gisaid.org) under the accession codes listed in Supplementary File 1. All codes for our analyses are available at https://github.com/AMC-LAEB/nl_sars-cov-2_genomic_epi_2022, copy archived at swh:1:rev:0a3c0a3e7d3959587e82e743162e28b45fa42dd7.

References

  1. Report
    1. Herrebrugh C
    2. Sluimer J
    3. Goderski G
    4. Benschop K
    5. Molenkamp R
    6. Meijer A
    7. Han W
    (2021)
    Report on research into the quality of SARS-CoV-2 antigen diagnostics in the Netherlands
    National Institute For Public Health and The Environment.

Decision letter

  1. Sarah E Cobey
    Reviewing Editor; University of Chicago, United States
  2. Miles P Davenport
    Senior Editor; University of New South Wales, Australia

Our editorial process produces two outputs: i) public reviews designed to be posted alongside the preprint for the benefit of readers; ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Regional importation and asymmetric within-country spread of SARS-CoV-2 variants of concern in the Netherlands" 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 Miles Davenport as the Senior Editor. The reviewers have opted to remain anonymous.

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:

1) Please streamline the narrative arc and logic of the paper. The reviewers questioned (including in the consultation session) what the estimation of variant growth rates adds to the analysis. It seems important to emphasize that the qualitatively estimated impacts of flight bans are if anything underestimated here, given other concurrent NPIs; please clarify which restrictions were in place. Please also see Reviewer 2's comment on the role of nightclubs.

2) There were technical concerns about the soundness of the results. It appears as though only one chain was used; several are traditional and safer. Biases might be introduced by the use of a strict vs. relaxed clock. A discrete phylogeographic model would be more appropriate. These improvements would lend greater confidence to the results. There should also be some treatment (at least in the Discussion) of the potential impact of sampling biases (inside and outside The Netherlands) on the results.

3) The reviewers found many of the figures confusing or unhelpful. Please revise them in light of the suggestions below.

Reviewer #1 (Recommendations for the authors):

1. I realize this is pretty foundational to phylogeographic methods (and welcome being referred to a textbook or classic paper I have missed), but how sensitive are the results – e.g., estimates of importation from ground travel/various areas---to sampling biases? For this question, I'm mostly thinking about deliberate subsampling biases (e.g., 2:1 Netherlands to non-Netherlands).

2. I'm curious about the potential magnitude of these effects due to biases in case ascertainment that might be poorly understood. In the U.S., there were many-fold differences in populations' access to testing and probability of being reported as a case, conditional on infection. These differences had strong socioeconomic and spatial associations. It seems deeply misleading to me to describe the sequencing as "random." It's not, unless case identification was random, which it wasn't. Maybe you could say something like "sequences were collected from a random sample of cases." But it would be good in the Methods and maybe the main text to describe exactly what is known about case ascertainment biases in The Netherlands over the study period, and the potential impacts of case ascertainment biases within The Netherlands and between countries on the estimated flows.

3. There's one mention of flight bans but otherwise "restrictions" is used. It would be useful to describe the restriction/ban on air travel precisely for people less familiar with the policies.

4. The github repo needs better documentation and structure to ensure reproducibility.

Reviewer #2 (Recommendations for the authors):

I have found the manuscript in need of several improvements:

1. The focus of the manuscript is somewhat diffuse and at times misguided, where for example the analysis of SARS-CoV-2 introductions into the Netherlands seems fine, others are questionable like using continuous instead of discrete phylogeography, and the analysis of VOC growth rate differences feeling a bit pointless, especially in light of the narrative the authors seem to push regarding the role of nightclubs in SARS-CoV-2 transmission. The authors should rethink the structure of the manuscript and what they want to say with it. If it's finding that travel restrictions do not achieve their goals – perfect! But adding analyses to beef up the manuscript afterward doesn't add much to the story.

2. There are a couple of technical issues that should be addressed. Firstly, all Bayesian MCMC analyses describe a single chain which is not standard. Ideally, at least two independent runs are expected for any MCMC analysis. Secondly, some of the figures don't seem like they convey much information, e.g. Figure 1C, Figure 3, and the bottom panel in Figure 4, and their use should therefore be reconsidered.

Overall, the manuscript has potential but requires substantial effort to improve, starting with its scope.

As stated in the public review part I think the authors should reframe, and refocus their manuscript. I think a study on the ineffectiveness of targeted flight cancellations is worthwhile but I don't see the added value (or good justification) of using continuous (instead of discrete) phylogeography or the logistic regressions.

I should also say that the narrative being pushed about nightclub reopenings is very flimsy and potentially confounded with the Netherlands having an age-tiered vaccination programme with age groups being given access around June 2021. This is not mentioned at all in the text and offers an alternative hypothesis for the distribution of cases across age groups during the reopening. I understand the authors are officially not claiming that there is strong evidence for it but when it comes up more than once in the text and is part of the discussion I personally would like it to be backed up with more evidence and all caveats addressed properly.

Line 76 "The four aforementioned VOCs also emerged in the Netherlands" – 'emerged' is often understood as 'originated', better to say 'arrived to' or 'dominated' here.

Figure 1A – Middle panel should say "Lineage proportions" for clarity, consider labelling each dominant area with the lineage that it represents since the legend in Figure 1C contains too many colours.

Figure 1A – The bottom panel should say "Mobility change (percent)" to make it intuitive.

Figure 1B – I have a personal preference for neatly-delineated colour maps and it's up to the authors if they want to change the limits between incidence and sequences so they're, for example, multiples of 20s and 500s, respectively.

Figure 1C – Not sure the tree is adding much here.

Line 143 "and continued to accumulate in frequencies" – change to 'frequency'.

Figure 2 – Worth shading European states with good sequencing programmes with a similar hue? Are those countries not flying passengers in?

Figure 3 – This figure is so messy perhaps it's not worth using? If the authors are dead-set on having some representation of migrations maybe there are better ways of doing it rather than a mess that says "there's a lot of movement"?

Figure 4 – Panel B is not labelled.

Geographic diffusion analyses indicate that more populous regions in the Netherlands play a large role in virus introduction and dissemination but surely this is entirely expected, given that such areas host transportation hubs for introductions and have more hosts to seed other areas.

Line 278 "locatthed" -> "located".

Paragraph 2 of the discussion contains passages that are more suitable for the Results section.

A very large dataset is analysed in BEAST so why did the authors decide to not use a relaxed clock and infer the rate naively? I find it very hard to believe that a dataset spanning close to a year would not have a temporal signal to inform the rate.

In all MCMC analyses, only one MCMC chain is mentioned. I'm sure the authors are aware that most studies are expected to run at least two independent ones, particularly with complex data like these.

Reviewer #3 (Recommendations for the authors):

1. Line 53: "Coronavirus-19 disease" should be "coronavirus disease 2019".

2. Line 88: samples were randomly selected for sequencing, not genomes.

3. Variants of concern (VOC), variants, and genotypes are used interchangeably. It would help to just pick of these and use them consistently.

4. I have a few suggestions regarding the visualization of Figure 1:

4.1. I would suggest using the Pangolin lineages as labels in panel C, and only list variants that are the focus of the study (all other lineages can be listed as "other").

4.2. The colors for 20D, 20I, 21A, 21I, and 21J are very similar (same in Figure 4). I would suggest using more distinctive colors to make it easier to see which color represents each of the variants.

4.3. It would be helpful to indicate the first detection of Α, Β, Γ, and Δ in Figure 1a.

5. Line 129. Did S Gene Target Failure (SGTF) influence the first detection of Α?

6. Figure 3. Panels B and D are hard to read, and may not be necessary to include. The advantage of only showing panels A and C is that their size can be increased.

7. Discussion. It would be helpful to start the discussion with a summary of the main findings. As raised before, I am not sure if the effect of international travel restrictions can really be distinguished from other interventions. Moreover, as Omicron was not part of the current study, I would refrain from speculating on the impact of those targeted flight restrictions as it is outside of the scope of the current study.

8. Limitations of the study are missing in the discussion.

9. Line 278: "locathed" should be "located".

10. Line 310. Methods lack a description of the used nucleic acid extraction method and diagnostic assay. Methods should also include a description of the used bioinformatics pipeline with specifics for how consensus genomes were generated (e.g. minimum depth and minimum frequency threshold to call consensus).

https://doi.org/10.7554/eLife.78770.sa1

Author response

Reviewer #1 (Recommendations for the authors):

1. I realize this is pretty foundational to phylogeographic methods (and welcome being referred to a textbook or classic paper I have missed), but how sensitive are the results – e.g., estimates of importation from ground travel/various areas---to sampling biases? For this question, I'm mostly thinking about deliberate subsampling biases (e.g., 2:1 Netherlands to non-Netherlands).

To clarify, we did not subsample at a 2:1 non-Netherlands to Netherlands sequence ratio for our global/Netherlands phylogenetic analyses. Instead, the number of sequences sampled per country per week is proportional to the number of reported cases that week. We ensured that we have at least one sequence from each country in each week. We aimed to sample around 100 sequences each week (the number may differ because of available sequences and the requirement to have at least one sequence from each country each week). As our sampling strategy may result in less non-Netherlands sequences than Netherlands sequences for some VOCs each week, we would then resample the non-Netherlands sequences again with a larger sample size, ensuring that we have at least twice as many non-Netherlands as Netherlands sequence that week.

Overall, we aim to be as comprehensive as possible to have all countries represented each week in our sampling strategy. We have now made this clearer in our methods:

Line 464: “As proportions of cases for each VOC lineage are unknown for most countries, we subsampled global sequences downloaded from GISAID (https://www.gisaid.org; dataset up to 6 October 2021) by the proportion of COVID-19 cases reported per week for each country using data from the Johns Hopkins University, Center for Systems Science and Engineering (CSSE) (http://github.com/CSSEGISandData/COVID-19). We aimed to sample 100 sequences each week. To ensure that countries that are underreporting cases (Gill et al., 2022) were included in our analyses, at least one representative sequence was included for each country with reported cases that week. The actual number of sequences sampled each week for each variant may differ because of sequence availability and the requirement to have at least one sequence from each country each week. We also subsampled Dutch sequences based on the weekly number of cases in different GGD regions as described above. As there may be certain weeks when our aforementioned strategy sampled less non-Netherlands sequences than Netherlands sequences, we would resample a larger number of non-Netherlands sequences accordingly such that there are at least twice as many non-Netherlands sequences as Netherlands sequences each week. The subsampling procedure yielded 6,365 (2,369), 1,531 (90), 1,274 (102) and 6,929 (1,035) Α, Β, Γ and Δ global (Dutch) sequences respectively.”

To check if our results are robust to sampling, we repeated our analyses with an independent sample bootstrap and found similar results as before:

Line 360: “We repeated our overseas phylogeographic analyses with an independent bootstrap subsample of sequence data (Figure 2 —figure supplements 2-3). There are differences in the estimated distributions of countries with nearest overseas neighboring sequences to Dutch subtrees (Figure 2 —figure supplement 1 and 3). However, our conclusions that there were multiple introductions from other European countries before and after travel restrictions were imposed on countries where VOCs first emerged remain evident in the bootstrap analysis (Figure 2 —figure supplement 2).”

2. I'm curious about the potential magnitude of these effects due to biases in case ascertainment that might be poorly understood. In the U.S., there were many-fold differences in populations' access to testing and probability of being reported as a case, conditional on infection. These differences had strong socioeconomic and spatial associations. It seems deeply misleading to me to describe the sequencing as "random." It's not, unless case identification was random, which it wasn't. Maybe you could say something like "sequences were collected from a random sample of cases." But it would be good in the Methods and maybe the main text to describe exactly what is known about case ascertainment biases in The Netherlands over the study period, and the potential impacts of case ascertainment biases within The Netherlands and between countries on the estimated flows.

We have now included further details in the Methods on the sources of the samples and how they were sampled for analyses:

Line 426: “These samples were collected through nationwide community testing programs including at test buses and health facilities (https://www.rivm.nl/coronavirus-covid-19/onderzoek/kiemsurveillance). Test-positive samples were randomly subsampled each week in each municipality to minimize case ascertainment bias. Only specimens with cycle threshold (Ct) values < 30 were selected for whole genome sequencing.”

3. There's one mention of flight bans but otherwise "restrictions" is used. It would be useful to describe the restriction/ban on air travel precisely for people less familiar with the policies.

We have now included more details on the travel restriction imposed on countries where VOCs first emerged:

Line 227: “To deter the introductions of novel VOCs into the Netherlands, travel restrictions were imposed on countries where the VOCs first emerged, including the United Kingdom between December 2020 and March 2021 due to the emergence of Α; South Africa and Brazil between January and June 2021 due to Β and Γ respectively; and India from April to June 2021 due to Δ. These travel restrictions include a ban on all incoming passenger flights except for those carrying cargo and medical personnel, on top of an entry ban for all non-European Union residents (Government of the Netherlands, 2021). On the other hand, travel within the Schengen Area of the European Union, which includes the Netherlands, remained possible during this period.”

4. The github repo needs better documentation and structure to ensure reproducibility.

We have now provided more detailed documentation in the GitHub repository to reproduce our results and figures.

Reviewer #2 (Recommendations for the authors):

I have found the manuscript in need of several improvements:

1. The focus of the manuscript is somewhat diffuse and at times misguided, where for example the analysis of SARS-CoV-2 introductions into the Netherlands seems fine, others are questionable like using continuous instead of discrete phylogeography, and the analysis of VOC growth rate differences feeling a bit pointless, especially in light of the narrative the authors seem to push regarding the role of nightclubs in SARS-CoV-2 transmission. The authors should rethink the structure of the manuscript and what they want to say with it. If it's finding that travel restrictions do not achieve their goals – perfect! But adding analyses to beef up the manuscript afterward doesn't add much to the story.

As stated in the Introduction, this work aims to identify where novel SARS-CoV-2 VOCs are likely introduced into the Netherlands and then characterize the spread of these VOCs in the country after their introduction. We selected different methods for our phylogeographic analyses based on the resolution of data available. We opted to perform discrete phylogeography to identify where likely overseas importation events were occurring from since only country/regional-level of geographic data was available to us. On the other hand, as postcodes of sampled Dutch individuals were available to us and since within-country spread can be assumed to occur mostly in a geographically homogenous way across the country, we employed continuous phylogeography to characterize the within-country spread of the novel variants. Continuous phylogeography methods have also been used to characterize within-country SARS-CoV-2 spread dynamics in several countries, including the United Kingdom (Kraemer et al., Science 2021; McCrone et al., medRxiv 2022), Rwanda (Butera et al., Nat Comm 2021) and Belgium (Bollen et al., Nat Sci Rep 2021).

We have removed the logistic regression analyses on VOC growth rate differences as we agree it adds little to the narrative.

2. There are a couple of technical issues that should be addressed. Firstly, all Bayesian MCMC analyses describe a single chain which is not standard. Ideally, at least two independent runs are expected for any MCMC analysis.

We had performed two independent Bayesian MCMC chains for our BEAST analyses. The maximum clade credibility trees were generated by combining the posterior distribution of trees from both chains. This is now explicitly stated in the Methods:

Line 592: “We performed two independent chains, each with 300 million MCMC generations for each variant analysis, sampling every 50,000 steps. The first 100 million steps were discarded as burn-in. Assessment of convergence (effective sample size > 200) for each chain was performed using Tracer v1.71 (Rambaut et al., 2018). The posterior distribution of trees from both chains were combined and resampled at every 100,000 steps to generate the maximum clade credibility tree.”

Secondly, some of the figures don't seem like they convey much information, e.g. Figure 1C, Figure 3, and the bottom panel in Figure 4, and their use should therefore be reconsidered.

We have removed Figure 1C and the bottom panel for Figure 4. We have revised and simplified Figure 3 to show the early within-country source-sink dynamics of Α and Δ during the first four weeks after their initial detection in the Netherlands.

Overall, the manuscript has potential but requires substantial effort to improve, starting with its scope.

As stated in the public review part I think the authors should reframe, and refocus their manuscript. I think a study on the ineffectiveness of targeted flight cancellations is worthwhile but I don't see the added value (or good justification) of using continuous (instead of discrete) phylogeography or the logistic regressions.

See response to comment 1 above.

I should also say that the narrative being pushed about nightclub reopenings is very flimsy and potentially confounded with the Netherlands having an age-tiered vaccination programme with age groups being given access around June 2021. This is not mentioned at all in the text and offers an alternative hypothesis for the distribution of cases across age groups during the reopening. I understand the authors are officially not claiming that there is strong evidence for it but when it comes up more than once in the text and is part of the discussion I personally would like it to be backed up with more evidence and all caveats addressed properly.

We agree and have removed all sections suggestive of nightlife reopening as the main driving force behind the spike in Δ transmissions in the Netherlands during July 2021.

Line 76 "The four aforementioned VOCs also emerged in the Netherlands" – 'emerged' is often understood as 'originated', better to say 'arrived to' or 'dominated' here.

We have now rephrased to: “The four aforementioned VOCs were also introduced into the Netherlands…”.

Figure 1A – Middle panel should say "Lineage proportions" for clarity, consider labelling each dominant area with the lineage that it represents since the legend in Figure 1C contains too many colours.

We have changed the y-label of the middle panel to “Lineage proportions”. We have also simplified

Figure 1A – The bottom panel should say "Mobility change (percent)" to make it intuitive.

The y-label is now “Mobility change (%)”.

Figure 1B – I have a personal preference for neatly-delineated colour maps and it's up to the authors if they want to change the limits between incidence and sequences so they're, for example, multiples of 20s and 500s, respectively.

We have now changed the bins to multiples of 20s or 1000s.

Figure 1C – Not sure the tree is adding much here.

As mentioned above, Figure 1C is now removed.

Line 143 "and continued to accumulate in frequencies" – change to 'frequency'.

Corrected.

Figure 2 – Worth shading European states with good sequencing programmes with a similar hue? Are those countries not flying passengers in?

We could not further delineate between European countries with different sequencing intensity in Figure 2 as the discrete phylogeographic analyses were performed at the continental level. However, we have now provided this stratification in the figure supplements 1 and 3 to Figure 2 for the nearest neighbouring overseas sequences to Dutch sequences. Flights were possible between most European countries during the study period albeit with varying vaccination and/or quarantine requirements.

Figure 3 – This figure is so messy perhaps it's not worth using? If the authors are dead-set on having some representation of migrations maybe there are better ways of doing it rather than a mess that says "there's a lot of movement"?

We have now revised Figure 3 to make it more informative to support our findings.

Figure 4 – Panel B is not labelled.

Figure 4B has been removed.

Geographic diffusion analyses indicate that more populous regions in the Netherlands play a large role in virus introduction and dissemination but surely this is entirely expected, given that such areas host transportation hubs for introductions and have more hosts to seed other areas.

We agree and that is why we performed the within-country continuous phylogeographic analyses to confirm our null expectations.

Line 278 "locatthed" -> "located".

Corrected.

Paragraph 2 of the discussion contains passages that are more suitable for the Results section.

We have now shifted those passages into the Results section, specifically:

Line 180: “As SARS-CoV-2 prevalence declined over time, average nationwide mobility also increased steadily which eventually came close to pre-pandemic levels in June 2021 (mean percentage change relative to pre-pandemic baseline = -5.0% (s.d. = 11.0%); Figure 1A).”

A very large dataset is analysed in BEAST so why did the authors decide to not use a relaxed clock and infer the rate naively? I find it very hard to believe that a dataset spanning close to a year would not have a temporal signal to inform the rate.

We redid our BEAST analyses using a relaxed lognormal clock.

Line 583: “We used HKY+G nucleotide substitution model and a Skygrid coalescent model (Hill and Baele, 2019) (each grid point denoting one week) with a relaxed lognormal molecular clock.”

In all MCMC analyses, only one MCMC chain is mentioned. I'm sure the authors are aware that most studies are expected to run at least two independent ones, particularly with complex data like these.

See response to comment 2 above.

Reviewer #3 (Recommendations for the authors):

1. Line 53: "Coronavirus-19 disease" should be "coronavirus disease 2019".

Corrected.

2. Line 88: samples were randomly selected for sequencing, not genomes.

Corrected.

3. Variants of concern (VOC), variants, and genotypes are used interchangeably. It would help to just pick of these and use them consistently.

We have now solely used “variants of concern (VOC)” throughout the text.

4. I have a few suggestions regarding the visualization of Figure 1:

4.1. I would suggest using the Pangolin lineages as labels in panel C, and only list variants that are the focus of the study (all other lineages can be listed as "other").

In response to comment 3 by reviewer 2, we have now removed Figure 1C.

4.2. The colors for 20D, 20I, 21A, 21I, and 21J are very similar (same in Figure 4). I would suggest using more distinctive colors to make it easier to see which color represents each of the variants.

We have now simplified the figure by distinguishing only the pre-Α and VOC lineages by different colours.

4.3. It would be helpful to indicate the first detection of Α, Β, Γ, and Δ in Figure 1a.

This is now included in Figure 1A.

5. Line 129. Did S Gene Target Failure (SGTF) influence the first detection of Α?

No.

6. Figure 3. Panels B and D are hard to read, and may not be necessary to include. The advantage of only showing panels A and C is that their size can be increased.

In response to comment 3 by reviewer 2, we have revised Figure 3 to make it more informative to support our findings.

7. Discussion. It would be helpful to start the discussion with a summary of the main findings. As raised before, I am not sure if the effect of international travel restrictions can really be distinguished from other interventions. Moreover, as Omicron was not part of the current study, I would refrain from speculating on the impact of those targeted flight restrictions as it is outside of the scope of the current study.

See response to comment 1 above.

8. Limitations of the study are missing in the discussion.

We have now expanded the Discussion on the limitations with our study:

Line 354: “Due to disparities in global genomic surveillance efforts and the lack of travel history information among the sampled Dutch individuals, we could not make more precise and accurate phylogeographic inferences on overseas introduction of VOCs into the Netherlands (Lemey et al., 2020). Furthermore, there were countries that are underrepresented or missing in our subsampling of non-Netherlands sequences. While we ensured that at least one sequence from each country with confirmed cases and genomic data was included in our analyses (see Methods), there may be overseas introductions that we could not account for. We repeated our overseas phylogeographic analyses with an independent bootstrap subsample of sequence data (Figure 2 —figure supplements 2-3). There are differences in the estimated distributions of countries with nearest overseas neighboring sequences to Dutch subtrees (Figure 2 —figure supplement 1 and 3). However, our conclusions that there were multiple introductions from other European countries before and after travel restrictions were imposed on countries where VOCs first emerged remain evident in the bootstrap analysis (Figure 2 —figure supplement 2).”

9. Line 278: "locathed" should be "located".

Corrected.

10. Line 310. Methods lack a description of the used nucleic acid extraction method and diagnostic assay. Methods should also include a description of the used bioinformatics pipeline with specifics for how consensus genomes were generated (e.g. minimum depth and minimum frequency threshold to call consensus).

We now provide the required information in the Methods:

Line 432: “As testing and samples were analyzed by 30-35 different laboratories across the country, different nucleic acid extraction methods were used (Herrebrugh et al., 2021). For samples analyzed by the laboratory of the Dutch National Institute for Public Health and Environment (38,260 samples; 96% of all samples analyzed), total nucleic acid was extracted using MagNApure 96 (MP96) with total nucleic acid kit small volume (Roche). RT-qPCR was performed on 5μl total nucleic acid using TaqMan Fast Virus 1-Step Master Mix (Thermo Fisher) on Roche LC480 II thermal cycler with SARS-like β coronavirus (Sarbeco) specific E-gene primers and probe and EAV as described previously (Corman et al., 2020; Scheltinga et al., 2005).

Amplicon-based SARS-CoV-2 sequencing for was performed using the Nanopore protocol “PCR tiling of COVID-19 virus (Version: PTC_9096_v109_revE_06FEB2020)” which is based on the ARTIC v3 amplicon sequencing protocol (Artic Network, 2020). Several modifications were made to the protocol as primer concentrations were increased from 0.125 to 1 pmol for the following amplicon primer pairs: in pool A amplicons 5, 9, 17, 23, 55, 67, 71, 91, 97 and in pool B amplicons 24, 26, 54, 64, 66, 70, 74, 86, 92, 98. Both libraries were generated using native barcode kits from Nanopore SQK-LSK109 (EXP-NBD104, EXP-NBD114 and EXP-NBD196) and sequencing was performed on a R9.4.1 flow cell multiplexing 2 up to 96 samples per sequence run. Consensus sequences (>50x depth-of-coverage) are generated using an in-house bioinformatic pipeline developed by the Dutch National Institute for Public Health and Environment (https://github.com/RIVM-bioinformatics/SARS2seq/).”

https://doi.org/10.7554/eLife.78770.sa2

Article and author information

Author details

  1. Alvin X Han

    Department of Medical Microbiology and Infection Prevention, Amsterdam University Medical Centers, Amsterdam, Netherlands
    Contribution
    Conceptualization, Formal analysis, Supervision, Validation, Investigation, Visualization, Methodology, Writing – original draft, Project administration, Writing – review and editing
    For correspondence
    x.han@amsterdamumc.nl
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6281-8498
  2. Eva Kozanli

    Department of Medical Microbiology and Infection Prevention, Amsterdam University Medical Centers, Amsterdam, Netherlands
    Contribution
    Formal analysis, Validation, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
  3. Jelle Koopsen

    Department of Medical Microbiology and Infection Prevention, Amsterdam University Medical Centers, Amsterdam, Netherlands
    Contribution
    Formal analysis, Validation, Investigation, Visualization, Methodology, Writing – original draft, Writing – review and editing
    Competing interests
    No competing interests declared
  4. Harry Vennema

    Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Bilthoven, Netherlands
    Contribution
    Resources, Data curation, Methodology, Writing – review and editing
    Competing interests
    No competing interests declared
  5. RIVM COVID-19 molecular epidemiology group

    Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Bilthoven, Netherlands
    Contribution
    Data curation, Formal analysis, Investigation, Methodology, Resources, Software, Validation, Writing – review and editing
    Competing interests
    No competing interests declared
    1. Lynn Aarts, Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Bilthoven, Netherlands
    2. Sanne Bos, Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Bilthoven, Netherlands
    3. Annemarie van den Brandt, Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Bilthoven, Netherlands
    4. Sharon van den Brink, Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Bilthoven, Netherlands
    5. Jeroen Cremer, Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Bilthoven, Netherlands
    6. Kim Freriks, Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Bilthoven, Netherlands
    7. Ryanne Jaarsma, Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Bilthoven, Netherlands
    8. Dennis Schmitz, Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Bilthoven, Netherlands
    9. Euníce Then, Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Bilthoven, Netherlands
    10. Bas van der Veer, Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Bilthoven, Netherlands
    11. Lisa Wijsman, Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Bilthoven, Netherlands
    12. Florian Zwagemaker, Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Bilthoven, Netherlands
  6. Karim Hajji

    Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Bilthoven, Netherlands
    Contribution
    Resources, Data curation, Software, Methodology
    Competing interests
    No competing interests declared
  7. Annelies Kroneman

    Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Bilthoven, Netherlands
    Contribution
    Resources, Data curation, Methodology
    Competing interests
    No competing interests declared
  8. Ivo van Walle

    Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Bilthoven, Netherlands
    Contribution
    Resources, Data curation, Methodology
    Competing interests
    No competing interests declared
  9. Don Klinkenberg

    Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Bilthoven, Netherlands
    Contribution
    Resources, Data curation, Methodology, Writing – review and editing
    Competing interests
    No competing interests declared
  10. Jacco Wallinga

    Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Bilthoven, Netherlands
    Contribution
    Resources, Data curation, Methodology, Writing – review and editing
    Competing interests
    No competing interests declared
  11. Colin A Russell

    Department of Medical Microbiology and Infection Prevention, Amsterdam University Medical Centers, Amsterdam, Netherlands
    Contribution
    Conceptualization, Supervision, Funding acquisition, Project administration, Writing – review and editing
    Competing interests
    No competing interests declared
  12. Dirk Eggink

    1. Department of Medical Microbiology and Infection Prevention, Amsterdam University Medical Centers, Amsterdam, Netherlands
    2. Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Bilthoven, Netherlands
    Contribution
    Conceptualization, Resources, Data curation, Supervision, Funding acquisition, Methodology, Project administration, Writing – review and editing
    Competing interests
    No competing interests declared
  13. Chantal Reusken

    Centre for Infectious Disease Control, National Institute for Public Health and the Environment, Bilthoven, Netherlands
    Contribution
    Conceptualization, Resources, Data curation, Supervision, Funding acquisition, Methodology, Project administration, Writing – review and editing
    For correspondence
    chantal.reusken@rivm.nl
    Competing interests
    No competing interests declared

Funding

European Research Council (ERC NaviFlu (no. 818353))

  • Alvin X Han

National Institutes of Health (5R01AI132362-04)

  • Colin A Russell

Nederlandse Organisatie voor Wetenschappelijk Onderzoek (09150182010027)

  • Colin A Russell

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank the administrators of the GISAID database for supporting rapid and transparent sharing of genomic data during the COVID-19 pandemic and all our colleagues sharing data on GISAID. A full list acknowledging the authors submitting genome sequence data used in this study can be found in Supplementary file 2. AXH and CAR were supported by ERC NaviFlu (no. 818353). C.A.R. was also supported by NIH R01 (5R01AI132362-04) and an NWO Vici Award (09150182010027).

Ethics

Human subjects: The Centre for Clinical Expertise at the National Institute for Public Health and the Environment (RIVM) assessed the research proposal following the specific conditions as stated in the law for medical research involving human subjects. The work described was exempted for further approval by the ethical research committee. Pathogen surveillance is a legal task of the RIVM and is carried out under the responsibility of the Dutch Minister of Health, Welfare and Sports. The Public Health Act (Wet Publieke Gezondheid) provides that RIVM may receive pseudonymized data for this task without informed consent. All necessary patient/participant consent has been obtained and the appropriate institutional forms have been archived, and any patient/participant/sample identifiers included were not known to anyone (e.g., hospital staff, patients or participants themselves) outside the research group so cannot be used to identify individuals.

Senior Editor

  1. Miles P Davenport, University of New South Wales, Australia

Reviewing Editor

  1. Sarah E Cobey, University of Chicago, United States

Publication history

  1. Received: March 18, 2022
  2. Preprint posted: March 22, 2022 (view preprint)
  3. Accepted: August 14, 2022
  4. Version of Record published: September 13, 2022 (version 1)

Copyright

© 2022, Han et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 271
    Page views
  • 58
    Downloads
  • 0
    Citations

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

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Alvin X Han
  2. Eva Kozanli
  3. Jelle Koopsen
  4. Harry Vennema
  5. RIVM COVID-19 molecular epidemiology group
  6. Karim Hajji
  7. Annelies Kroneman
  8. Ivo van Walle
  9. Don Klinkenberg
  10. Jacco Wallinga
  11. Colin A Russell
  12. Dirk Eggink
  13. Chantal Reusken
(2022)
Regional importation and asymmetric within-country spread of SARS-CoV-2 variants of concern in the Netherlands
eLife 11:e78770.
https://doi.org/10.7554/eLife.78770
  1. Further reading

Further reading

    1. Epidemiology and Global Health
    Carlos A Prete Jr, Lewis Fletcher Buss ... Ester C Sabino
    Research Article

    Background: The COVID-19 situation in Brazil is complex due to large differences in the shape and size of regional epidemics. Understanding these patterns is crucial to understand future outbreaks of SARS-CoV-2 or other respiratory pathogens in the country.

    Methods: We tested 97,950 blood donation samples for IgG antibodies from March 2020 to March 2021 in eight of Brazil’s most populous cities. Residential postal codes were used to obtain representative samples. Weekly age- and sex- specific seroprevalence was estimated by correcting the crude seroprevalence by test sensitivity, specificity and antibody waning.

    Results: The inferred attack rate of SARS-CoV-2 in December 2020, before the Gamma VOC was dominant, ranged from 19.3% (95% CrI 17.5% - 21.2%) in Curitiba to 75.0% (95% CrI 70.8% - 80.3%) in Manaus. Seroprevalence was consistently smaller in women and donors older than 55 years. The age-specific infection fatality rate (IFR) differed between cities and consistently increased with age. The infection hospitalisation rate (IHR) increased significantly during the Gamma-dominated second wave in Manaus, suggesting increased morbidity of the Gamma VOC compared to previous variants circulating in Manaus. The higher disease penetrance associated with the health system's collapse increased the overall IFR by a minimum factor of 2.91 (95% CrI 2.43 - 3.53).

    Conclusions: These results highlight the utility of blood donor serosurveillance to track epidemic maturity and demonstrate demographic and spatial heterogeneity in SARS-CoV-2 spread.

    Funding: This work was supported by Itaú Unibanco 'Todos pela Saude' program; FAPESP (grants 18/14389-0, 2019/21585-0); Wellcome Trust and Royal Society Sir Henry Dale Fellowship 204311/Z/16/Z; the Gates Foundation (INV- 034540 and INV-034652); REDS-IV-P (grant HHSN268201100007I); the UK Medical Research Council (MR/S0195/1, MR/V038109/1); CAPES; CNPq (304714/2018-6); Fundação Faculdade de Medicina; Programa Inova Fiocruz-CE/Funcap - Edital 01/2020 Number: FIO-0167-00065.01.00/20 SPU Nº06531047/2020; JBS - Fazer o bem faz bem.

    1. Epidemiology and Global Health
    2. Medicine
    David Robert Grimes
    Tools and Resources

    There is increasing awareness throughout biomedical science that many results do not withstand the trials of repeat investigation. The growing abundance of medical literature has only increased the urgent need for tools to gauge the robustness and trustworthiness of published science. Dichotomous outcome designs are vital in randomized clinical trials, cohort studies, and observational data for ascertaining differences between experimental and control arms. It has however been shown with tools like the fragility index (FI) that many ostensibly impactful results fail to materialise when even small numbers of patients or subjects in either the control or experimental arms are recoded from event to non-event. Critics of this metric counter that there is no objective means to determine a meaningful FI. As currently used, FI is not multi-dimensional and is computationally expensive. In this work a conceptually similar geometrical approach is introduced, the ellipse of insignificance (EOI). This method yields precise deterministic values for the degree of manipulation or miscoding that can be tolerated simultaneously in both control and experimental arms, allowing for the derivation of objective measures of experimental robustness. More than this, the tool is intimately connected with sensitivity and specificity of the event / non-event tests, and is readily combined with knowledge of test parameters to reject unsound results. The method is outlined here, with illustrative clinical examples.