Abstract
The COVID-19 pandemic has posed an ongoing global crisis, but how the virus spread across the world remains poorly understood. This is of vital importance for informing current and future pandemic response strategies.
We performed two independent analyses, travel network-based epidemiological modelling and Bayesian phylogeographic inference, to investigate the intercontinental spread of COVID-19.
Both approaches revealed two distinct phases of COVID-19 spread by the end of March 2020. In the first phase, COVID-19 largely circulated in China during mid-to-late January 2020 and was interrupted by containment measures in China. In the second and predominant phase extending from late February to mid-March, unrestricted movements between countries outside of China facilitated intercontinental spread, with Europe as a major source. Phylogenetic analyses also revealed that the dominant strains circulating in the USA were introduced from Europe. However, stringent restrictions on international travel across the world since late March have substantially reduced intercontinental transmission.
Our analyses highlight that heterogeneities in international travel have shaped the spatiotemporal characteristics of the pandemic. Unrestricted travel caused a large number of COVID-19 exportations from Europe to other continents between late February and mid-March, which facilitated the COVID-19 pandemic. Targeted restrictions on international travel from countries with widespread community transmission, together with improved capacity in testing, genetic sequencing and contact tracing, can inform timely strategies for mitigating and containing ongoing and future waves of COVID-19 pandemic.
Introduction
The coronavirus disease 2019 (COVID-19), caused by the highly contagious causative agent severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), was first identified in Wuhan, Hubei province, in late 2019.1,2 The World Health Organization declared COVID-19 a pandemic on 11 March 2020. As of 1 October 2020, there have been more than 33 million confirmed COVID-19 cases with 1.01 million fatalities worldwide, affecting more than 200 countries, territories or areas.3 Non-pharmaceutical interventions are currently the only viable strategies available to contain and suppress COVID-19 spread,4–6 especially for the travel-associated spread between and within countries.7–13 For example, Wuhan’s lockdown on 23 January 2020 delayed the growth and limited the size of the COVID-19 outbreak in China.7–9 These travel restrictions also had an effect in reducing international transmission of SARS-CoV-2, with a 77% reduction in cases exported from mainland China to other countries.10,11 Accordingly, epidemics triggered by international importations in multiple countries have been suppressed by substantial travel and social distancing interventions.12,13 However, widespread community transmission of SARS-CoV-2 was still reported in some regions after easing these measures, with an increase in confirmed cases or subnational localized resurgences in the USA, India, Japan and some European countries since late June.14,15
Studies have attempted to quantify the spread risk of COVID-19 using epidemiological models within a country and at regional scales,5,10,16–18 and accumulated divergences among genomic data of SARS-CoV-2 have also been used to trace the viral spread. For instance, it was found that most of the SARS-CoV-2 strains circulating in the NY area and northern CA were introduced from Europe via multiple independent importation events.19,20 Similarly, the great majority of the SARS-CoV-2 strains circulating in UK was derived from continental Europe in March via inbound travellers.21 The effects of international travel and border control measures on the global spread of SARS-CoV-2 were also preliminarily explored.10,22 However, questions remain as to the spatiotemporal characterizations and transmission routes of the COVID-19 pandemic and which factors contributed to the seeding of the virus and the emergence of outbreaks. Answers to these questions are of vital importance in formulating effective measures to contain the ongoing COVID-19 pandemic and for future infectious disease outbreak planning.
Using a comprehensive dataset integrating epidemiological, travel, intervention and genetic data, here we conduct epidemiological modelling and Bayesian phylogeographic inference to: (i) understand the changing patterns of international movements under COVID-19 interventions across the world, (ii) measure the transmission dynamics of intercontinental COVID-19, (iii) distinguish the main source of international transmission that facilitated the COVID-19 pandemic and further (iv) reveal the source of the dominant strains circulating in the USA. The findings of our study can be used to inform public health response efforts for ongoing and future waves across the world.
Methods
Simulating the COVID-19 spread using an epidemiological model
Using domestic and international population mobility data and a susceptible-exposed-infectious-removed (SEIR) framework,6 we built a global travel network-based stochastic metapopulation transmission model to simulate COVID-19 transmission across 221 countries, territories or areas globally from 1 December 2019 to 31 March 2020 (the R code of the SEIR model can be found at: https://github.com/wpgp/BEARmod). To initially parameterize the model, country-specific reproduction numbers [median: 2.4, interquartile range (IQR): 2.0–2.8], before the implementation of physical distancing interventions, were estimated from daily case counts reported by each country, adjusted for reporting delays.23 A publicly available dataset of government COVID-19 countermeasures was obtained to define the timings of various interventions implemented in each country.24 We also used initial epidemiological information (i.e. the incubation period and delays from symptom onset to report) estimated from case data during the early stages of the outbreak in Wuhan.25
Additionally, two mobile phone-derived publicly available anonymized population mobility datasets covering 133 countries, territories or areas in 202026,27 were used in our simulations to account for the impact of travel and physical distancing interventions implemented to mitigate the COVID-19 pandemic across space and time. To simulate the international spread of COVID-19 through population movement, we parameterized the model using global air traffic data, from 1 December 2019 to 31 March 2020, obtained from the Official Aviation Guide (https://www.oag.com/). As we simulated the COVID-19 spread by day, we used the proportion of daily volume of seats on scheduled flights over the total number of flight seats in each month and route to disaggregate air passenger counts from monthly to daily level between countries. Details of the model, parameters and data sources and their collation are provided in the Supplementary Materials and Methods section.
Bayesian phylogeographic inference of SARS-CoV-2
All available complete or near-complete genomic nucleotide sequences of SARS-CoV-2 were collected from GISAID (https://www.gisaid.org), NCBI GenBank (https://www.ncbi.nlm.nih.gov), NMDC (http://nmdc.cn/coronavirus) and NGDC (https://bigd.big.ac.cn/ncov/) databases on 1 April 2020. The latest collection date of the sequences was 26 March 2020. We excluded duplicated sequences, genomes with incomplete collection date and low-quality sequences. We kept one representative of any known epidemiologically linked case clusters (http://virological.org/t/phylodynamic-analysis-176-genomes-6-mar-2020/356). Further, we removed sequences isolated from Africa and South America due to a small number of available sequences as of 1 April 2020. We pooled the remaining genomes into five geographic regions, including China, the Asian countries minus China, North America, Europe and Oceania. There were 14 jackknife resampled datasets generated for the phylogeographic analyses by randomly sampling at most 10 isolates per location per day. Each down-sampled dataset comprised a total of 1030 sequences, including China (n = 233), the Asian countries without China (n = 118), North America (n = 311), Europe (n = 296) and Oceania (n = 72). Virus sequences were aligned with default parameters in MAFFT v7.28
To investigate the global phylodynamics of SARS-CoV-2, phylogeographic inferences were performed using the asymmetric discrete trait analysis in BEAST v1.10.4 (Supplementary Materials and Methods section).29 We employed the general time reversible (GTR) model to describe nucleotide substitution and to account for rate heterogeneity across sites by combining the discrete gamma distribution with four rate categories and invariable sites. A strict molecular clock was utilized to model sequence evolution. Two tree priors, exponential growth and Bayesian skyline, were employed to describe the virus population dynamics. By evaluating the marginal likelihood and evolutionary rate estimates, the exponential growth model was selected in subsequent analyses. Also, we set up Markov jump count parameters to obtain the frequency of transitions between location traits along phylogenetic branches over time (https://beast.community/markov_jumps_rewards). Analyses with the same parameter configuration were performed on the 14 subsampled datasets to assess the robustness of phylogeographic inferences. We performed each analysis with 100 million iterations, later removing an appropriate burn-in (10–40% of samples) to get an adequate effective sample size (basically ESS ≥100). The maximum clade credibility (MCC) tree was calculated from the posterior tree space in TreeAnnotator v1.10.4.
Phylogenetic and bioinformatic analyses of worldwide SARS-CoV-2
To investigate which lineage of SARS-CoV-2 contributed to the epidemic in different locations, we performed a maximum likelihood analysis of the global SARS-CoV-2 strains (n = 5482 and updated on 14 April 2020 in GISAID) using RAxML v8.2.930 with 100 bootstrap replicates and the general time reversible (GTR) model. We classified global SARS-CoV-2 into two lineages based on two highly linked single nucleotide polymorphisms (SNPs): type S (8782T and 28144C) and type L (8782C and 28 144T).31,32
Based on the maximum likelihood analysis and the specific nucleotide mutations highlighted on the Nextstrain website (https://nextstrain.org/), we proposed a simple classification frame that was able to capture the major clades of global SARS-CoV-2. In detail, the lineage L was further classified into L1 (241C, 3037C, 14408C, 23403A) and L2 (241T, 3037T, 14408T, 23403G) (Supplementary Materials and Methods section). The lineage L2 was further classified into L2.1 (28881–28883AAC) and L2.2 (25563T), which has a special sub-lineage L2.2.1 (1059T). Analogously, type S was further divided into S1 (18060C, 17858 and 17747) and S2 (18060T, 17858G and 17747T). The numbering of nucleotide sites is referred to the NCBI reference sequence of SARS-CoV-2, NC_045512. Further, SARS-CoV-2 strains isolated from the USA by 20 May 2020 were classified into different lineages. The distributions of different virus lineage within each state of the USA were summarized.
Results
Changing patterns of international travel
Compared with the levels of air travel between 1 December 2019 and 22 January 2020, international outbound travel from China following the lockdown of the country rapidly declined to a median of 16% (IQR: 9–55%). However, international air travel from countries outside of China did not significantly change in February, and reductions for European and Middle Eastern countries began to appear in early March. During the week of 25–31 March 2020, international air travel across the world reduced to a median of 37% (IQR: 18–63%) of the levels seen before the pandemic (Supplementary Figures 1 and 2).
Two distinct phases of intercontinental COVID-19 spread
Two transmission phases of international exportations of COVID-19 were identified by both the epidemiological model and Bayesian phylogeographic inference (Figure 1). During the first phase in mid-to-late January, 519 cases [95% credible interval (CI): 369–682] were estimated to have been exported from China to other countries, with most of them (84%) occurring before 1 February 2020. However, 3 weeks later, more intercontinental exportations of COVID-19 occurred in the second phase, from late February to mid-March (Figure 1A and Supplementary Figure 3). We estimated that there were 5563 cases (95% CI: 2714–8627) exported from Asian countries, except China, to other continents by 31 March. Europe was estimated to have the highest number of exported cases (17 087, 95% CI: 11 342–23 136) to other continents, with > 90% of exported cases occurring in March. North America and Oceania were estimated to have exported 3462 cases (95% CI: 2169–4817) and 1428 cases (95% CI: 662–2231) to other continents, respectively. The decline of this transmission phase corresponded with the international travel reductions implemented in mid-to-late March (Supplementary Figures 1 and 2). Additionally, our estimated epicurves were consistent with the patterns seen for reported curves as of 31 March 2020 (P < 0.001, R2 = 0.96), and high correlations were also found between the reported data and the estimated incidence by country (P < 0.001, R2 = 0.70) and estimated imported cases from China (P < 0.001, R2 = 0.65), respectively (Supplementary Figures 4 and 5).
Dispersal dynamics of global COVID-19 cases and SARS-CoV-2 virus across regions. (A) The estimated numbers of COVID-19 cases exported from different regions by 31 March 2020. (B) Migration dynamics of SARS-CoV-2 virus from the different regions by 26 March 2020. The grey vertical dotted lines from left to right indicate the date of Wuhan shutdown in China, the date of COVID-19 declared as a Public Health Emergency of International Concern (PHEIC) and the date of COVID-19 declared as a pandemic by WHO, respectively. Colours in the plot represent the different regions defined in the present study.
We summarized the Markov jumps between defined locations in the phylogenies of the posterior tree space to describe the intercontinental transmissions of SARS-CoV-2, and the inferred phylodynamics of SARS-CoV-2 from 14 subsampling datasets are consistent and robust. Consistent with results from the SEIR modelling (Figure 1A), two peaks of global SARS-CoV-2 transmission were also identified by Bayesian phylogeographic inference (Figure 1B). The first peak of virus dissemination occurred during mid-to-late January 2020, but shortly after the Wuhan shutdown, viral transmission from China reduced dramatically. The second peak, which was higher than the first, occurred around early to mid-March and decreased rapidly after 15 March. The Bayesian inference also highlighted that China was the major contributor to the first spread peak, while multiple regions contributed to the second spread peak, with Europe being the dominant source.
The second phase facilitating the COVID-19 pandemic
We further studied the viral importation events for major continental regions. The travel network-based SEIR modelling revealed two distinct peaks of COVID-19 importation events in most continents by 31 March 2020 (Figure 2). We estimated that there were 450 cases (95% CI: 338–568) exported from China into other Asian countries in the first phase and 11015 cases (95% CI: 7008–15 235) from other continents into Asia in the second phase. Europe had 28 estimated cases (95% CI: 9–50) imported from China in the first phase and 5664 cases (95% CI: 2759–8832) from other continents in the second phase. Similarly, most of the imported cases in North America were from other continents (5559, 95% CI: 3447–7834) in the second phase, with only 23 cases (95% CI: 14–31) imported from China. Meanwhile, Africa, South America and Oceania were estimated to have had 3316 (95% CI: 1846–4948), 2077 (95% CI: 1299–2909) and 1114 (95% CI: 677–1584) cases, respectively, imported from other continents in the second phase, with few cases exported from China into Africa (5 cases, 95% CI: 1–12), South America (1 case, 95% CI: 0–2) and Oceania (13 cases, 95% CI: 7–20). However, travel restrictions implemented in mid-to-late March have significantly curbed the international spread of COVID-19 since the end of March.
The estimated numbers of COVID-19 cases imported into different continents. This estimation ends at March 31, 2020. The median and IQR of estimates are provided. The three grey dotted vertical lines are the same as those in Figure 1.
The Bayesian inference analyses also found that viral exportations from China peaked during mid-to-late January and dropped sharply after the shutdown of Wuhan (Figure 3A). In contrast, exportations from Europe rapidly increased from mid-February and peaked in early March. During this period, the virus migration frequency from Europe to North America and Oceania was estimated to have increased until the WHO declared COVID-19 a pandemic on 11 March (Figure 3B). Additionally, virus diffusions from Europe to North America were likely to have already occurred before travel restrictions on European citizens were implemented by the USA on 13 March 2020. However, since 26 March 2020, shortly after the second transmission peak, the USA has been the country with the highest reported number of COVID-19 cases.
Migration events of SARS-CoV-2 from China and Europe to other continents. (A) Estimated Markov jumps of SARS-CoV-2 from China to other regions by 26 March 2020. (B) Estimated Markov jumps of SARS-CoV-2 from Europe to other regions by 26 March 2020. The mean and 95% high density intervals of estimates are provided. A Markov jump represents the frequency of transitions between locations along phylogenetic branches over time and can be regarded as one migration event of SARS-CoV-2. The three grey dotted vertical lines are the same as those in Figure 1.
Origins of the dominant strains circulating in the USA
To further understand the impact of international travel on the introduction of SARS-CoV-2 in the USA, we further simulated importation events into the USA. Our epidemiological SEIR model revealed that most of the cases imported into the USA (as of 31 March 2020) likely originated from Europe (2673 cases, 95% CI: 1951–3401), other North American countries (1324 cases, 95% CI: 924–1730) and Oceania (192 cases, 95% CI: 89–297) (Figure 4A). However, only a limited number of cases (18, 95% CI: 12–23) were imported from China, with 685 cases (95% CI: 336–1040) coming from other Asian countries (Figure 4A).
Temporal dynamics of the imported COVID-19 cases and the circulating SARS-CoV-2 lineages in the USA. (A) Estimated COVID-19 case numbers imported into the USA from different sources (see legend) by 31 March 2020. Cases imported from China are highlighted in the inset box. The three grey dotted vertical lines are the same as those in Figure 1. (B) A simple classification framework for SARS-CoV-2 based on the maximum likelihood analysis and specific nucleotide substitutions along the SARS-CoV-2 genome. (C) Spatial distribution of the reported COVID-19 cases and the SARS-CoV-2 lineages within mainland USA. The shading from light to dark blue indicates the number of cumulative COVID-19 cases reported in each state by 20 May 2020. Names of the top 10 states with the highest cumulative COVID-19 cases are underlined. The size of each circle represents the number of SARS-CoV-2 genome sequences in each state available from GISAID by 20 May 2020. Colours in the circles represent proportions of different virus lineages
Clearly, several SARS-CoV-2 variants of both S and L genotypes have been co-circulating in the USA (Figure 4C). The clade L2.2.1, characterized by 241T, 3037T, 14408T, 23404G, 25563T and 1059T, was dominant (~51.07% of the US strains available in GISAID as of 20 May 2020) in most of the States of the USA, including at least 7 out of the top 10 states that have reported the most COVID-19 cases as of 20 May, particularly in NY (78.84%), NJ (84.00%) and MI (77.37%) (Figure 4C and Supplementary Table 1). In 13 out of the 14 MCC trees, L2.2.1 was predicted to have originated from Europe in late February (Supplementary Figure 6). Similarly, recent research also revealed that most of the SARS-CoV-2 strains circulating in the NY area and northern CA were introduced from Europe.20,33 Apart from L2.2.1, there were also independent introductions from Europe to North America in other clades (Supplementary Figure 6). Although an S genotype variant, S2 (characterized by 8782 T, 28144C and 18060T), was dominant (54.94%) in the WA state and was predicted to have originated from China, it accounted for a minority of the confirmed SARS-CoV-2 cases in the USA (Supplementary Figures 7 and 8 and Supplementary Table 1).34 The dominant SARS-CoV-2 strains contributing to the community transmission in the USA originated from Europe and were introduced before the implementation of travel restrictions to European citizens, followed by extensive domestic transmission.
Discussion
Using comprehensive and integrated datasets on population mobility, interventions, case reporting and virus genomes, we have conducted both epidemiological modelling and Bayesian phylogeographic inferences to uncover the intercontinental dynamics of COVID-19 transmission and virus evolution. The findings from the two independent methodologies were consistent and both revealed two phases of COVID-19 transmission across the globe. In the first phase, the virus mainly circulated and propagated in China, and our modelling only found a small number of COVID-19 cases transmitted from China to other localities. This first phase was interrupted by strict containment measures implemented proactively across China and reduced outbound travel.5,25
However, in the second intercontinental dispersal phase, the seeding of viruses was exaggerated by population movements between countries outside of China from late February through mid-March, which finally led to the global pandemic. In particular, our study identified Europe as an epicentre of this phase through phylogenetic analyses, which is consistent with other studies that revealed the dominant strains circulating in the USA and the UK before travel restrictions were imported from continental Europe where community transmission was widespread.19,21,35,36 The USA, as the most-affected country, had substantial airline network connections to Europe but imposed relatively late travel restrictions on European travellers. Our phylogenetic analyses revealed the dominant strains in 7 out of 10 states with the highest case numbers were imported from Europe in the second phase rather than from China in the first phase. This was in accordance with the results from epidemiological modelling and further highlighted the role of international travel from Europe in SARS-CoV-2 spread into the USA and across the globe.
More importantly, apart from Europe, several other regions outside of China, including North America and some Asian countries, were also estimated to have contributed to the international spread of SARS-CoV-2. However, the reductions in outbound travel from China since late January, as well as strict restrictions on international travel across the world since late March, have effectively reduced international spread of the virus.6–7,10,22,25 However, many countries where COVID-19 had been contained or eliminated are now facing increasing risks of international importations after relaxing travel and social distancing interventions. For example, more than 2000 internationally imported COVID-19 cases have been documented in China as of 25 July 2020.37 Although strict testing, quarantine and contact tracing measures have been implemented, some cases have caused new outbreaks in different regions of China, including Harbin, Shulan, Beijing, Dalian and Urumqi in May–July 2020, with resurgences also reported in the USA, India, Japan and several European countries.3,37–40 The relaxation of travel restrictions and social distancing, particularly in the summer holiday period, and the reopening of schools might have contributed to the resurgence and ongoing outbreaks of COVID-19 in Europe. The potential impact of restoring domestic and international mobility on the resurgence and subsequent international spread of COVID-19 warrants further study. Further, the coordination between countries and regions in relaxing interventions and lifting international travel restrictions can greatly improve the likelihood of containing and reducing COVID-19 spread post-lockdown.4
Our findings should be considered in the context of several assumptions and data limitations. First, we only estimated the spread risk of COVID-19 via air travel here due to data availability. With the integration of land and sea transport models, the international travel networks are more complex and our results may underestimate the intercontinental risk of COVID-19 spread via travellers. Second, the accuracy of our modelling relies on accurate estimates of epidemiological parameters partially derived from reported case data, the quality of which might be constrained by case definitions, diagnosis and surveillance capacity and other factors across countries/regions.41–43 Third, mobile phone-derived data for parameterizing travel and physical distancing interventions in our model may not be representative of the population in each country due to variations in coverage and the spatiotemporal distribution of users in populations. Fourth, other factors and interventions, such as hand washing and wearing facemasks, may also contribute to suppress COVID-19 spread or migration across space and time,44–47 but our simulations did not specify their contributions to international transmission. Finally, our Bayesian phylogeographic inferences only included SARS-CoV-2 genomes available in GISAID by 26 March 2020, and the potential sampling bias among different continents may result in underestimates of the peak of the second transmission phase and cause biases in the estimates.
However, our analyses highlight how heterogeneities in intercontinental travel have facilitated the intercontinental seeding of SARS-CoV-2 and shaped the spatiotemporal characteristics of the ongoing COVID-19 pandemic. Unrestricted movements between continents/countries after easing stringent lockdown measures would likely trigger a new wave of COVID-19 spread across the world. Our findings improve our understanding of early transmission dynamics of COVID-19 across continents and can help to tailor public health response strategies accordingly. Specifically targeted international travel interventions should be adapted to different phases of the pandemic for countries and their corresponding travellers,48 together with timely testing, genetic sequencing and contact tracing49 for COVID-19 infections in travellers, which are also needed to monitor the ongoing pandemic and mitigate resurgences post-lockdown.
Authors’ contribution
Conceptualization: W.S., Y.B. and S.L.; Data curation: J.Y., J.L. and S.L.; Formal analysis: J.Y., J.L., S.L. and P.W.; Visualization: J.Y., J.L. and S.L.; Writing: J.Y., S.L., Y.B., W.S., C.W.R., A.C., N.W.R., A.J.T., J.R.F. and W.X.; Original draft: J.Y., S.L., Y.B., W.S. and A.J.T. All authors read and approved the manuscript.
Acknowledgements
We thank the important work of SARS-CoV-2 genome data producers globally contributing sequence data to the GISAID, NCBI GenBank, NMDC and NGDC databases. We also acknowledge Google and Baidu for sharing publicly available population mobility data.
Conflict of interest
The authors have declared no conflicts of interest.
Funding
This work was supported by the Key Research and Development Project of Shandong Province (2020SFXGFY01, 2020SFXGFY08); National Key Research and Development Programme of China (2020YFC0840800, 2016YFE0205800); Strategic Priority Research Programme of the Chinese Academy of Sciences (XDB29010102, XDA19090118); National Natural Science Foundation of China (32041010, 81773498); National Science and Technology Major Project of China (2016ZX10004222-009); Academic Promotion Programme of Shandong First Medical University (2019QL006); Bill & Melinda Gates Foundation (OPP1134076); European Union Horizon 2020 (MOOD 874850); National Major Project for Control and Prevention of Infectious Disease in China (2017ZX10104001, 2018ZX10101004); National Natural Science Foundation of China Outstanding Young Scholars (31822055 to Y.B.); Youth Innovation Promotion Association of CAS (2017122 to Y.B.); Bill & Melinda Gates Foundation (OPP1170969 to N.W.R.; OPP1106427, OPP1032350, OPP1134076 and OPP1094793 to A.J.T.); Clinton Health Access Initiative, UK Department for International Development and Wellcome Trust (106866/Z/15/Z, 204613/Z/16/Z to A.J.T.); Young Taishan Scholars Program of Shandong Province of China (tsqn20161046 to W.X.); Academic Promotion Programme of Shandong First Medical University (2019RC010 to W.X.); Taishan Scholars Programme of Shandong Province (ts201511056 to W.S.); and China Postdoctoral Science Foundation (2020T130123ZX to J.Y.).
References
Author notes
Jing Yang, Juan Li and Shengjie Lai contributed equally to this work
