Editor's evaluation: Regional importation and asymmetric within-country spread of SARS-CoV-2 variants of concern in the Netherlands
Year: 2022
Type: peer-review
Abstract: Article Figures and data Abstract Editor's evaluation Introduction Methods Results Discussion Data availability References Decision letter Author response Article and author information Metrics 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 Decision letter Reviews on Sciety eLife's review process 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. Figure 1 Download asset Open asset 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 over time, average nationwide mobility also increased which eventually to pre-pandemic levels in June 2021 percentage change relative to pre-pandemic baseline Figure 1A). SARS-CoV-2 was at its then with reported positive cases that week. 1 week after however, weekly cases 50,000 on weeks 26 and 27 2021). most infections attributed to the novel VOC the Alpha VOC as the dominant lineage within the 3 weeks as over of randomly genomes were as Delta VOCs by the number of weekly reported positive cases by age the relative proportions in case positive remained the study period except for weeks 26 and 27 where the increase in cases was largely attributed to between 15 and (Figure 1A). of the transmission among then was linked to nightlife et al., 2022). In response, the government nightlife establishments down again on 2021 (week numbers after but remained at over new cases per week for the rest of the study period. The Delta VOC had in displaced Alpha by then with over of randomly genomes sampled from August 2021 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 including the United Kingdom between December 2020 and March 2021 due to the emergence of South and between January and June 2021 due to Beta and Gamma and India from April to June 2021 due to These travel restrictions include a on all except for and on of an for all of the Netherlands, 2021). On the other hand, travel within of the in which the Netherlands, remained largely during this period. To identify likely where and when VOCs were introduced into the country, we subsampled a representative set of Dutch and overseas sequences collected over the same time period. We then reconstructed maximum (ML) phylogenies and used these fixed trees to perform discrete analyses using a Bayesian to likely overseas introductions at the continental level. This was by Dutch sequences with ancestral that were attributed to an overseas (Figure All four VOCs were introduced into the Netherlands prior to the targeted flight restrictions that were imposed on countries where these VOCs first emerged. Figure 2 with 3 all Download asset Open asset Likely overseas introduction of VOC lineages into the Netherlands at the continental level. For each VOC a time-scaled maximum phylogeny using the Dutch and their overseas neighboring sequences was analyses were performed to the likely continental region of ancestral or with ancestral nodes attributed to an overseas but Dutch sequences are area the when a targeted flight was imposed on the country where the VOC lineage first (i.e. (A) Alpha, United (B) Beta, South Gamma, Importantly, countries where travel restrictions were in we estimated multiple likely introduction from other foreign countries into the Netherlands for all four VOCs Beta, Gamma, in global sequencing efforts et al., 2021), the random surveillance strategy used in local sample and genetic among SARS-CoV-2 genomes used to ancestral we are to and the number of introductions attributed to different geographical However, many of the estimated regions for these ancestral were in Europe of all estimated overseas introduction Beta, Gamma, Figure these European introductions to during the targeted travel period. the phylogenetic ancestral to the aforementioned we found that many of these overseas neighboring were detected in and where between the Netherlands remained as well as other countries U.K., where targeted travel restrictions were set in place (Figure There was also period in time in which these VOCs were introduced into the Netherlands introductions likely occurred during the period when these VOCs were also within the 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 downsampled sets of Alpha and Delta sequence data (Figure For the first four weeks the detection of both VOCs within the country, introductions and phylogenetic branch were in the more populous regions of the country, including North and South Holland, and a of early dominant this to regions of these GGD regions occurred as well but are less However, as local infections were seeded in these bidirectional in phylogenetic between different regions the country. These bidirectional continued to increase as of the VOC over time, a lockdown in the case of Alpha (Figure In for the Delta VOC, we observed a spike in spread the week of June 2021, with estimated increase in total phylogenetic branch by 6 2021 (week This in likely to the case numbers observed between weeks 25 and 27 June 2021). Figure 3 Download asset Open asset dynamics the Alpha and Delta variants of concern within the Netherlands during the first four weeks after their respective (A) Alpha 2 and December (B) Delta 20 April and 2021). Figure Download asset Open asset number of phylogenetic branch and rate of Alpha and Delta VOCs in the Netherlands. line the number of branch from early dominant regions including North and South Holland, and line number of branch from of these early dominant Discussion international travel restrictions are in the Netherlands is to importation risks of novel SARS-CoV-2 variants from its regional due to within the European As this regional is not to the Netherlands and has been reported in other European countries as well (Lemey et al., 2021; et al., 2021; et al., 2021). Importantly, regional introductions of novel lineages new waves of infections in Europe (Lemey et al., 2021). to September 2020, the dominant variant lineages (i.e. 20 A and 20E that circulated the Netherlands were seeded by from its European (Hodcroft et al., 2021; et al., 2021; Munnink et al., 2020). In the introduction of SARS-CoV-2 in the Netherlands in February 2020 were attributed to where the sustained European SARS-CoV-2 transmission was seeded (Oude Munnink et al., 2020; et al., 2020). Here, we that all four VOC lineages detected in the country up to August 2021 also from its European Importantly, regional importation risks the period these VOCs circulated the country and with where targeted flight restrictions were imposed on countries where these VOCs first emerged. The emergence of the VOC in in 2021 et al., 2022) again led to targeted flight restrictions by countries in the including the Netherlands which was a Delta infection However, the VOC was detected in samples collected one week before the imposed flight and did not to become the dominant VOC circulating in the Netherlands by the end of 2021 that travel restrictions are useful restrictions from most countries provided that local is in the first place et al., 2021). to in global genomic surveillance efforts and the of travel information among the sampled Dutch we not more and phylogeographic on overseas introduction of VOCs into the Netherlands (Lemey et al., 2020). there were countries that are or in our subsampling of non-Netherlands sequences. we that at least one sequence from each country with cases and genomic data was included in our analyses there may be overseas introductions that we not We our overseas phylogeographic analyses with an independent of sequence data (Figure There are in the estimated of countries with overseas neighboring sequences to Dutch (Figure 1 and However, our that there were multiple introductions from other European countries before and after travel restrictions were imposed on countries where VOCs first in the (Figure We also found that early introductions of VOCs, the Alpha and Delta VOCs, are more likely found in populous regions of the
Source: eLife
Author: Sarah E Cobey
Institution: University of Chicago
Language: English
Subfield: Modeling and Simulation
Field: Mathematics
Domain: Physical Sciences
Open Access status: gold
Estimated APC paid: $2,500
OpenAlex ID: W4297805154