The origin and early spread of SARS-CoV-2 in Europe

Introduction

In response to the pandemic potential of the SARS-CoV-2 virus, many nations closed their borders in order to curb the virus’ spread (1). These closures continue to incur high economic and social costs. To weigh the relative costs and benefits of border closures, it will be important to understand the efficacy of these policies. At the early stages of an outbreak, border closures can delay a pathogen’s arrival, thereby giving countries additional time to prepare (2). However, the success of this strategy depends on timely implementation and a good knowledge of where the pathogen is already circulating. To evaluate the efficacy of border closures in limiting the spread of SARS-CoV-2, it is important to reconstruct the timeline of the early international spread of the virus, before such policies were implemented.

In this analysis, we aim to estimate the early patterns of SARS-CoV-2 transmission into and across Europe. We also address the more specific question of where the predominant SARS-CoV-2 lineage circulating in Europe originated. We hope that by addressing these questions we can inform further analysis of the efficacy of border closures as a strategy to combat SARS-CoV-2.

The SARS-CoV-2 virus was identified as the cause of an epidemic in Wuhan, China in late 2019 (3). The epidemic in Wuhan was reported to the WHO on 31 Dec. 2019 and within one month, SARS-CoV-2 was confirmed to have spread to 19 additional countries (4). By the end of February 2020, the virus was detected in all WHO regions (5). Currently, several lineages of the SARS-CoV-2 virus are circulating across the globe. The intermixing of these lineages in different countries and regions suggests that the virus has been transmitted across borders many times (6).

Here we focus on estimating the early introductions of SARS-CoV-2 into Europe and the virus’ migration across European borders. Through national surveillance efforts, the first COVID-19 cases in Europe were detected in France on 24 Jan. 2020 and in Germany on 28 Jan. 2020 (7, 8). Of the 47 cases detected in Europe by 21 Feb. 2020, 14 were infected in China, 14 were linked to the initial cases in Germany, 7 were linked to the initial cases in France, and 12 were of unknown origin (7). In addition to the unknown sources of transmission, some early introductions may not have been detected. This is especially probable given that a significant proportion of infected individuals are likely to be asymptomatic (9). In summary, it is difficult to draw firm conclusions about the source, number, and timing of SARS-CoV-2 introductions into Europe based on confirmed case data alone.

Viral genomes are an important secondary source of information on outbreak dynamics. If viruses acquire mutations on the same timescale as an outbreak, these mutations can provide information about past transmission events. Phylodynamic methods couple a model of viral evolution describing the mutational process to an epidemiological model describing the transmission process. By fitting the combined model to viral genomes sampled from a cohort of infected individuals, we can infer the evolutionary and epidemiological model parameters. Here we fit a phylodynamic model with geographic structure to SARS-CoV-2 genomes from Hubei, China and several European countries before the first borders were closed in these regions. We co-infer the transmission tree linking these sequences, the geographic location of ancestral lineages, migration rates of infected individuals between regions, the effective reproductive number, and the proportion of no-longer infectious cases sequenced in each region.

In addition to these inferences, we specifically focus on estimating the geographic origin of the predominant SARS-CoV-2 lineage in Europe. This lineage is defined by a characteristic amino acid substitution at position 314 in the ORF1b gene from proline to leucine and was provisionally named the “A2a” lineage by the Nextstrain team. In the more dynamic, tree-based nomenclature suggested by (10), this lineage corresponds to the “B.1” lineage described as “a large lineage that roughly corresponds to the large outbreak in Italy, and has since seeded many different countries” (11). As of Apr. 1, 2020, two-thirds of the SARS-CoV-2 sequences collected in Europe belonged to this lineage and just 10% of sequences from the lineage were collected outside Europe (data from (12), lineages assigned using (13)). Here, we use the name A2a to refer to the group of SARS-CoV-2 viruses defined by the ORF1b:P314L mutation.

Where the A2a lineage originated remains unclear. Its characteristic ORF1b mutation was found in some of the earliest confirmed COVID-19 cases in Italy, Switzerland, Germany, Finland, Mexico, and Brazil in late February (14, 15). Intriguingly, a late-January sample from a German case linked to business travel from Singapore shares a mutation in the S gene with the A2a lineage, but does not have the lineage-defining ORF1b mutation. This German sample is part of a smaller “A2” clade that is basal to the larger clade of A2a sequences (6). As a result, it was hypothesized that a German transmission cluster may have seeded the larger European outbreak (14–16). However, it was quickly pointed out that incomplete and biased sampling must be taken into account before this hypothesis can be rigorously addressed (14, 17, 18).

Phylodynamic models with geographic structure aim to account for such biases. Firstly, parameter estimates are generated by integrating over a distribution of potential phylogenies, which acknowledges that we cannot reconstruct the true transmission tree with certainty. Secondly, sampling parameters are allowed to differ between regions, which acknowledges that testing and sequencing resources vary across regions. Here, we fit a phylodynamic model with geographic structure to full-length SARS-CoV-2 genomes to (i) estimate the early patterns of SARS-CoV-2 spread into and across Europe, (ii) weigh genomic evidence for competing hypotheses about the geographic origin of the predominant A2a lineage in Europe, (iii) report on the epidemiological parameters, and (iv) compare the rate of new cases arising from within-region transmission versus migration during the early epidemic.

Materials and Methods

Results

Inference results: A2a lineage origin

The maximum clade credibility tree in Figure 1 summarizes the posterior sample of transmission trees linking analyzed sequences. The A2a lineage sequences form a clear clade with posterior support of 1. The most recent common ancestor of the analyzed A2a sequences is estimated to be in Italy with 87% posterior support. In contrast, the location of the most recent common ancestor between this clade and the basal, Singapore-linked German sequence is less certain. This ancestor is inferred to have been in either Germany (40% posterior support), Hubei (38%), or Italy (20%). We find very little support for this ancestor having been in France or another European country (2%).

• Download figure
• Open in new tab

Figure 1.

Maximum clade credibility tree. The clade of A2a sequences analyzed is highlighted with dashed branches. The values above the branches are the posterior clade probabilities and the pale red bars show the 95% highest posterior density interval for node ages. The pie charts at nodes show posterior support for the ancestor being located in each deme (note that we assumed the root of the tree was in Hubei with probability 1). The deme for each tip is the deme in which the sequence was collected, irrespective of travel history. Tips are annotated with GISAID accession identifier.

Comparing rates of migration and within-region transmission

Figure 2 compares the rate at which we estimate new cases to arise in each region from migration versus from within-region transmission. The estimated rates of new cases from migration and within-region transmission are represented here as point estimates 5 days before the date of case confirmation, which assumes a 5-day delay between infection and onward transmission or migration. Beginning with the first day on which we have case data from Hubei, we estimate a substantial risk of infected individuals migrating from Hubei into European regions. Throughout late January to mid-February 2020, cases were sporadically detected in each European region, each of which is associated with a risk of subsequent within-region transmission. Sustained within-region transmission is first evident in Italy in mid-February. Shortly thereafter, sustained within-region transmission occurred in other European countries, in France, and in Germany. By 8 Mar. 2020, the estimated rate of occurrence of new cases from within-region transmission is within the estimated bounds on the rate of new cases from migration for each region considered (Figure S12A). We obtain the same qualitative result in our sensitivity analysis using a very different prior on the migration rate (Figure S12B). We note that the rates in Figure 2 are underestimates of the rates of new cases arising due to migration or transmission due to the underreporting in the confirmed case data. However, assuming that the amount of underreporting is comparable across regions, we can indeed compare the rates.

• Download figure
• Open in new tab

Figure 2.

Estimated rate of new cases arising from migration compared with the estimated rate of new cases arising from within-region transmission. For each day, we multiplied the number of newly confirmed cases in each source region by the posterior sample of migration rates from source to sink. The median of these rates is show in the “Migration” row. We also multiplied the number of newly confirmed cases in each sink region by the posterior sample of transmission rates for the region. The median of these rates is shown in the “Within-region transmission” row. Dates are lagged 5 days to account for a 5-day delay between infection and migration or onward transmission. Case data comes from (23).

Discussion

We inferred the early spread of the SARS-CoV-2 virus into and across Europe as well as the geographic origin of the predominant A2a lineage spreading in Europe. To do this, we applied a previously published phylodynamic model to analyze publicly available viral genome sequences from the epidemic origin in Hubei, China and from the earliest detected and largest European outbreaks before 8 Mar. 2020. After performing Bayesian inference, we (i) report on inferred patterns of SARS-CoV-2 spread into and across Europe, (ii) compare posterior support for several hypotheses on the origin of the A2a lineage, (iii) report on epidemiological parameters, and (iv) compare the timeline of new cases resulting from migration versus within-region transmission in Europe before borders were closed.

Genome sequence data indicates that prior to 8 Mar. 2020, SARS-CoV-2 was introduced from Hubei province into France, Germany, Italy and other European countries at least 2-4 times each (Table 1). These estimates, which are based solely on genome sequence data, provide a complementary account of introduction events compared to line-list data (26). The introduction events we report here are inferred to have occurred along the transmission tree specific to the analyzed sequence set and are not attributable to individual cases. In comparison, line-list data (7, 26) attributes introduction events to individual cases but cannot reconstruct previous, unobserved introductions. Since we analyze only a fraction of all cases, we expect our estimates to be a minimum bound on the true number of introductions.

Ideally, we want to go beyond counting migration events amongst the analyzed sequences and investigate general dynamics. To do this, we would interpret inferred migration rates as representing more general patterns of SARS-CoV-2 spread. However, the sequence data was only informative for inferring some of these rates (Figure S1). In regions with few lineages circulating during the period considered, there is little signal for the amount of outward migration. We observe information about the per-individual migration rate from Italy to other European countries (Figure S1). However, we do not find evidence for preferential migration paths from Hubei into different European regions or from each European region to the others, although we cannot exclude this possibility.

We estimate that the A2a viruses spreading in Europe by 8 Mar. 2020 had a common ancestor in Italy sometime between mid-January and mid-February 2020 (Figure 1). In contrast, Nextstrain places this ancestor in the U.K. with 100% confidence (27). However, the Nextstrain result may be an artefact of disproportionately high sequencing effort in the U.K. since biased sampling violates the assumptions of the “mugration” method employed (28). We additionally report that the A2a lineage was most likely carried from Hubei to Italy or from Hubei to Italy via Germany. Both transmission routes have almost equal posterior support under our model assumptions (Figure 1). Since we only consider a few geographic regions, these migration routes are not necessarily comprehensive. Rather than reconstructing a complete transmission chain, we compare model support for different A2a transmission routes amongst the analyzed demes and report two equally plausible routes.

Although it is not the main focus of our analysis, we also report on epidemiological parameters of the early outbreaks considered. Estimates for the reproductive number fall roughly within the range of previous estimates (29), though we mention a particular caveat with respect to the reproductive number in Hubei below. Unsurprisingly, prevalence estimates in early March generally exceed confirmed case counts by a factor of 1-3 (Figure S4). Our inferences of epidemiological parameters do not challenge the idea that the early reproductive number in different outbreaks is difficult to estimate precisely, but not hugely variable, and that there is substantial under-reporting in line-list data (30).

Finally, we estimated the rate of new cases arising from migration compared with the rate of new cases arising from within-region transmission in the regions analyzed. The magnitudes of these rates are quite uncertain due to uncertainty in the inferred migration and transmission rates (Figure S11) and under-reporting in case counts, which we implicitly assume to be constant in time and between demes. However, the temporal trends suggested by these data are still compelling and robust towards different prior assumptions. We see that under sustained risk of case migration from Hubei, isolated cases were confirmed throughout Europe beginning in late January 2020 but did not immediately cause large outbreaks. Shortly after the first evidence of sustained within-region transmission in Italy, outbreaks in the rest of Europe also took hold (Figure 2).

Our results based on the multi-type birth-death model take into account phylogenetic uncertainty and sampling biases between demes, which are two major concerns in genomic analyses of SARS-CoV-2 (18). Indeed, wide confidence intervals around internal nodes in the maximum clade credibility tree and low clade support near the tips (Figure 1) indicate a high degree of phylogenetic uncertainty. Therefore, it is important that the parameter estimates we report result from integrating over a distribution of potential phylogenies with different geographic locations assigned to ancestral lineages. In comparison, some initial studies that estimated international SARS-CoV-2 spread constructed a median-joining network instead of a phylogeny to account for this uncertainty (16, 31). In this approach, identical sequences are collapsed to single nodes and edges represent mutational differences. This disregards information from relative sampling times and means that ancestor-descendent relationships are highly dependent on the choice of the network root (32, 33). Unaccounted-for sampling biases in these analyses may also yield spurious results for the geographic origin of lineages (34, 35). Our analysis, which relies on a mechanistic model of migration and between-deme sampling differences, should be robust to such biases.

Despite the advantages of the multi-type birth-death model just mentioned, there are also several unique caveats to consider. The birth-death model assumes uniform-at-random sampling from the total infected population in each deme. However, particularly in the early stages of outbreaks, infected individuals were identified by health ministries via contact tracing (7). Non-random sampling may be one possible explanation for why we infer markedly different transmission rates in China when analyzing cases from within Hubei (as in this analysis) as opposed to cases exposed in Hubei but sequenced elsewhere (as in our previous analysis (36)). Furthermore, the multi-type birth-death model assumes that parameters are constant through time and homogenous within demes. As a result, our inferences based on province-, country-, and continent-level demes are only coarse approximations of the true, heterogeneous epidemic dynamics occurring at a local level. In particular, we do not account for a reduction in airline traffic between Europe and China beginning in late January, before borders were closed (37). Due to these limitations, we focus on estimating and interpreting particular events along the transmission tree of the analyzed sequences (e.g. Table 3, Figure 1) and advise caution when interpreting inferred migration rates (Figure S1).

We expect that our results will be useful in parameterizing more specialized models aimed to understand the efficacy of border closures as a means to fight pandemic disease. So far, such analyses have primarily used line-list data and information on travel networks to estimate SARS-CoV-2 migration patterns (38–40). Here we present independent estimates of migration patterns based on genome sequence data. By combining case count data and our estimates for migration and transmission rates, we provide a timeline of early SARS-CoV-2 introduction and spread before border closures were implemented. Despite migration risk from Hubei being on the same order of magnitude as later migration risk from Italy, we only observe sustained outbreaks in other European regions after the onset of sustained within-region transmission in Italy. Finally, before the first border closures in Europe, we estimate the risk of new cases arising from within-region transmission to be within the estimated range for the risk of new migration cases.