Establishment & lineage dynamics of the SARS-CoV-2 epidemic in the UK

Genomic data

All SARS-CoV-2 genomes available on GISAID (14) on 23 June 2020 were downloaded and combined with all SARS-CoV-2 genomes sequenced by the COG-UK consortium (15) by 26 June 2020 (available at https://www.cogconsortium.uk/data/). The pipeline used to collect and process raw SARS-CoV-2 sequence data and sample-associated metadata across the national COG-UK network is described in (37). We removed sequences that were from duplicate or environmental samples, those without exact collection dates, and those with large clusters of substitutions or large indels. Each genome sequence was aligned to the reference (Wuhan-Hu-1, GenBank: MN908947.3) using minimap v2.17 (38) and the resulting SAM alignment was converted to a FASTA alignment, with the 5’ and 3’ UTRs of each genome masked by Ns. Insertions relative to the reference were discarded and site 11,083 (site position relative to MN908947), which is globally homoplasic, was also masked. Genomes that contained >5% Ns after mapping and those with a genetic distance to WH04 (GISAID: EPI_ISL_406801) more than 4 standard deviations from the epi-week mean genetic distance to WH04 were discarded. The final dataset consisted of 50,887 genomes sampled between 24 December 2019 and 22 June 2020, of which 26,181 (∼51%) were from the UK (see Fig. 1A).

Geographical metadata

Administrative level 2 (admin2) metadata for the sampling location of UK virus genome sequences in the dataset (roughly equivalent to counties in the UK) required cleaning in order to be mapped to official admin2 regions, as found in the Global Administrative Database (GADM, https://gadm.org).

Some sampling locations in the metadata could not be unambiguously mapped to a known location (e.g. “City Centre”), while others were for locations in overseas territories (e.g. Falklands and Gibraltar). Yet other genome sequences had uninformative spatial records (e.g. Yorkshire or Wales), or no admin2 level data at all. For these (3431 of 26,181) the admin2 region was not mapped. We carried out a simple one-to-one mapping where possible, which included correcting spelling mistakes and alternative entries for the same county (e.g. Durham versus County Durham). Locations recorded at a higher spatial resolution were mapped to the corresponding admin2 region (e.g. Solihull was mapped to Birmingham). Where the recorded locations were larger than the admin2 regions (e.g. “West Midlands”), and most of the sequences in the area were from this larger conglomeration as opposed to its higherresolution components, these admin2 regions were combined. When creating the map figures, we also merged some city authorities with no reported sequences with their surrounding county, on the assumption that the larger county was used to represent the location of city samples (e.g. for Leicester and Leicestershire). Finally, genome sequences from Northern Ireland reported locations as historical counties, rather than the official admin2 designations, and so these historical counties were used instead. The cleaning code is provided on the GitHub repository (https://github.com/COG-UK/UK-lineage-dynamics-analysis).

Phylogenetic analysis and molecular clock dating

We developed a new Bayesian molecular clock phylogenetic analysis pipeline in order to reconstruct a posterior set of time-scaled phylogenetic trees for our exceptionally large virus genome dataset. Using the standard Bayesian approach it is currently impractical to estimate time-scaled trees directly from genome sequence data for more than a few thousand sequences. Therefore, we employed a number of extensions to make the analysis tractable.

First, we divided the full genome sequence dataset (n = 50,887) into five smaller datasets. Genomes were assigned SARS-CoV-2 lineages according to the nomenclature defined in (39) using Pangolin (40; github.com/covlineages/pangolin). Each lineage (and its sublineages) represents a monophyletic clade in the global SARS-CoV-2 phylogeny and can thus be analysed independently. For each lineage in A (n = 3591), B (n = 8821), B.1 (n = 22,861), B.1.1 (n = 15,616), we estimated an approximately maximum-likelihood tree using the Jukes-Cantor model in FastTree v2.1.10 (41), then collapsed branch lengths shorter than 5×10⁻⁶ substitutions per site, which corresponded to distances smaller than one substitution across the whole virus genome, and likely result from nucleotide ambiguity codes in the genome sequences. By pruning out a large monophyletic clade the maximumlikelihood tree for B.1 was further divided into two trees, B.1.pruned (n = 12,275) and B.1.X (n = 10,586).

Prior to analysing the full dataset, an initial analysis was performed on a subset of genomes to obtain estimates of the molecular clock rate and of the TMRCA of each large-scale phylogenetic tree defined above. The full dataset was subsampled as evenly as possible across epi-weeks and countries with a slight enrichment for samples immediately descended from five large polytomies in the global phylogeny. For each of these nodes, we always included the five oldest genomes, the most recent genome sequence and five other immediate descendants that were randomly chosen. The remaining genomes were sampled by allocating an even number of sequences per epi-week while maintaining a dataset size of <1,000 genomes. For each epi-week, genomes were sampled evenly by country until either its allocation was exhausted or there were no remaining genomes available. This subsampled dataset was analysed in BEAST 1.10 (42) using a GTR+G+F substitution model, with a strict molecular clock model using a non-informative continuous-time Markov chain (CTMC) prior (43) and a Skygrid coalescent tree prior (44) with 40 grid points, roughly corresponding to weeks between 1 October 2019 and 2 July 2020. In the analysis, monophyly constraints were used to ensure that the clades corresponding to the large-scale phylogenetic trees identified in the previous step were monophyletic. We combined four independent Markov Chain Monte Carlo (MCMC) chains that were each run for 40 million steps, discarding the first 4 million steps of each chain as burn-in and resampling states every 4000 steps. Convergence was assessed using Tracer (45).

Next, we applied a commonly used approach, recently implemented in BEAST 1.10, to convert branches of the large-scale phylogenetic trees from units of substitutions per site to time. This model takes the place of the nucleotide substitution model in a traditional Bayesian molecular clock dating analysis. Briefly, each branch of a maximum-likelihood tree is first scaled to represent the number of substitutions that occurred along that branch. Polytomies are resolved by inserting branches of length 0 substitutions. The likelihood of a branch b_i of length s_i substitutions is defined by a Poisson distribution with mean t_im where t_i is the length of the branch in years and m is the clock rate. The log-likelihood of the whole tree is then the sum of the log-likelihoods of each branch, which represents a fixed, strict-clock model and follows a commonly implemented approach for scaling phylogenies into time-calibrated trees (e.g. 46-48).

Each large-scale phylogenetic tree was analysed under a strict clock model, with the clock rate fixed to the median estimate from the preliminary analysis (7.5×10⁻⁴ substitutions/site/year) and a Laplace root-height prior with mean equal to the median TMRCA estimate of the corresponding subtree in the preliminary analysis and scale equal to the average distance from the median. Trees were sampled using MCMC under the model described above with a Skygrid coalescent tree prior (44) using the same grid-points as in the preliminary analysis. A randomly resolved time-calibrated tree estimated in TreeTime (49) was used as the starting tree. To maintain a mapping between the topology in the estimated time-calibrated tree and the input genetic distance tree, we constrained the topologies such that any tree-move that broke a clade present in the input tree was rejected. The resulting MCMC chain, therefore, only samples different polytomy resolutions and branch durations. This approach allowed us to incorporate uncertainty in the polytomy resolutions and branch durations into our molecular clock analysis.

We ran between 8 and 24 chains for 60 to 100 million MCMC steps for each large-scale phylogeny. Upon completion, we discarded 15 million states as burn-in from each chain. Chains that did not converge or pass the burn-in in less than 15 million states were re-run. Chains were combined and resampled every 100,000 states using custom R-scripts, leaving between 6808 and 17,020 posterior samples of each large-scale phylogenetic tree. Convergence was assessed using Tracer (45) and the R-package coda (50).

Identifying transmission lineages

We define a “UK transmission lineage” as two or more UK infection cases that (i) descend from a shared, single importation of the virus into the UK from elsewhere, (ii) are the result of subsequent local transmission within the UK, and (iii) were present in our virus genome sequence dataset. This concept is illustrated in Figure S1 and is distinct from a transmission cluster, an epidemiological term commonly referring to a group of cases that occur close to each other in space and time (e.g. in a hospital or care home). Therefore, a large UK transmission lineage may comprise many different individual transmission clusters.

[It is important to note that the “UK transmission lineage” definition employed here is distinct from the lineage/phylotype designations used by other parts of the COG-UK consortium and that are displayed at https://microreact.org/project/cogconsortium. Those latter designations (which have the format “UK…”) are defined on the basis of shared sets of mutations, rather than shared descent from an inferred single introduction event.]

We can identify UK transmission lineages in the time-calibrated trees estimated in the previous step as clades of two or more genomes sampled in the UK. The TMRCA of all genome sequences in a UK transmission lineage represents the earliest transmission event in the lineage revealed by the data; however, it does not necessarily represent the first transmission event in the lineage as a whole, nor does it represent the importation date (i.e. the arrival date of the index patient in the UK). The relationship between the TMRCA of a UK transmission lineage in our dataset and the importation date is illustrated in Figure S2. Specifically, if the transmission lineage is well-sampled, then the TMRCA represents the date of the first transmission event in the lineage (TMRCA A in Fig. S2, UK transmission lineage 2 in Fig. S1). However, if the transmission lineage is sparsely sampled then the TMRCA may represent a later transmission event (TMRCA B in Fig. S2, UK transmission lineage 1 in Fig. S1). The “importation date” of each UK transmission lineage is the date that an infected inbound traveller entered the UK.

We used a two-state asymmetric discrete trait analysis (DTA) model (51) implemented in BEAST 1.10 (42) to infer ancestral node locations (UK, non-UK) on empirical distributions of 500 time-calibrated trees sampled from each of the posterior tree distributions estimated above. Additionally, we used a robust counting approach (52) to estimate the expected number of location state transitions into and out of the UK. For each large-scale subtree, we combined 2 independent chains, each run for 5 million MCMC steps and sampled every 4500 states. The first 10% of each run was discarded as burn-in, resulting in 2000 trees with estimates of the ancestral location for each internal node. Finally, TreeAnnotator 1.10 was used to generate maximum clade credibility (MCC) trees for each subtree, where each internal node is assigned a posterior probability of representing a transmission event in the UK.

Transmission lineages were identified by first labelling each node in the MCC trees as UK or non-UK and then initiating a depth-first search from each UK genome in the MCC trees. All nodes with a median age after 23 January 2020 and posterior probability >0.5 of the ancestral location being located in the UK were labelled as UK nodes. The depth-first search is continued until a non-UK node is encountered or there are no nodes left to explore. At the end of the depth-first search, all nodes visited by the search are added to the same (arbitrarily named) UK transmission lineage. If only one tip is visited, the UK genome at the tip is marked as a singleton. This procedure is repeated iteratively until every UK genome in the tree has been assigned a transmission lineage or marked as a singleton. The same procedure was repeated on each of the 2000 posterior trees, for each subtree, from the DTA analyses described above to examine statistical uncertainty in the number, size and duration of UK transmission lineages and their TMRCA distribution.

Our methodology is likely to underestimate the true number of transmission lineages and singletons. Since only a small fraction of UK infections have been sequenced (Fig. 1A), many lineages will have gone undetected. Furthermore, the power to detect a transmission lineage in our sparsely-sampled dataset is dependent on its size (i.e. the frequency of a lineage being sampled from a small random sample of infections), making it more likely for larger lineages to be detected. The low sampling fraction means that some singletons detected in our dataset likely belong to observed and unobserved UK transmission lineages. Nonetheless, the true number of singletons (importations not resulting in onward transmission) is likely to be significantly more than our estimate, because their small size makes them difficult to detect with a low sampling fraction. Finally, under-sampling of genomes from other countries could result in mistaken aggregation of separate importations, reducing the number of detected lineages. This mistaken aggregation will result in larger, older lineages being estimated. This was the motivation for placing an age limit on UK nodes in the tree. We chose 23 January 2020 as the oldest possible date for a transmission event in the UK as this represents the date that the first patient who tested positive for SARS-CoV-2 in the UK entered the country (53) (tested positive on 30 January 2020). Although older importations into the UK could in theory be possible, if they had resulted in large autochthonous outbreaks, we would have observed this in both epidemiological and genomic data.

We estimate a median of 2968 (95% HPD 2829-3103) non-UK to UK state transitions and an additional 1468 (95% HPD 1362-1566) UK to non-UK state transitions (Fig. S3, Table S1) using the robust counting approach (52). The former slightly exceeds the sum of transmission lineages and singletons as identified on the MCC tree (=2918) and across the 2000 posterior trees (median=2829, 95% HPD=2773-3048; Table S1). This result is expected, since multiple location state changes along long branches contribute to the total number of state transitions, but do not add to the total number of UK transmission lineages or singletons. The largest number of location state transitions occur on the B.1.1 phylogeny, with the fewest occurring on lineage A, which are the largest and smallest of the subtrees, respectively. Proportional to the number of tips, fewer state changes are inferred on the two B.1 phylogenies than other subtrees, while the number of UK to non-UK transitions on the B phylogeny exceeds that inferred on other lineages. We caution that UK to non-UK transitions are likely to be underestimated because of under-sampling in other countries and differences in the proportion of infections sequenced between countries.

The transmission lineage size distribution from the MCC trees falls within the HPD interval taken across the 2000 posterior trees (Fig. 1B). Although the sizes of the largest transmission lineages vary substantially across posterior trees, the cumulative size distributions are similar across all trees (Fig. 1B, inset). Similarly, the transmission lineage duration distribution on the MCC trees falls within the variation of the HPD interval taken across the 2000 posterior trees (Fig. S8).

We used the Jaccard index to compare the classification of UK genome sequences into transmission lineages and singletons between posterior trees and the MCC trees. Figure S13A shows the mean, median and 95% HPD interval of the Jaccard index for each posterior tree compared to the 1999 other posterior trees, across all subtrees. While most Jaccard indices are between 0.7 and 0.8, there is a noteworthy minority of trees with mean Jaccard indices <0.6 (n=100). Comparing the 2000 posterior trees to the classification on the MCC tree (Fig. S13B), results in a similar distribution of Jaccard indices, with most indices between 0.7 and 0.8 and minorities below 0.6 and above 0.8 (n=68, n=170 respectively).

We undertook a similar analysis of the sensitivity to phylogenetic uncertainty of the distribution of UK transmission lineage TMRCAs. We computed the median and 95% HPD interval of the number of transmission lineage TMRCAs on each date across the 2000 sampled posterior trees. Figure S14 shows that the TMRCA distribution computed from the MCC trees falls within the comparatively narrow HPD limits, and oscillates around the median estimate for each date.

UK epidemiological data

The number of reported COVID-19 cases in the UK, by specimen date, were downloaded from https://coronavirus.data.gov.uk/cases (date accessed: 1 September 2020). The number of reported COVID-19 cases for each Upper Tier Local Authority (UTLA) in England, Local Health Board (LHB) in Wales and regional NHS Board in Scotland, by specimen date, were downloaded from https://coronavirus.data.gov.uk/downloads/csv/coronavirus-cases_latest.csv, http://www2.nphs.wales.nhs.uk:8080/CommunitySurveillanceDocs.nsf (file: “Rapid COVID-19 surveillance data.xlsx”) and https://github.com/DataScienceScotland/COVID-19-Management-Information (file: “COVID19 - Daily Management Information - Scottish Health Boards - Cumulative cases.csv”), respectively (date accessed: 15 October 2020).

To enable comparison of case and sequence data, locations used to report case data were combined to correspond to those used for sequence data and vice-versa (see the Geographical metadata section). Northern Ireland was not included due to inconsistencies between the locations used for case and sequence data reporting that could not be easily resolved.

Global deaths due to COVID-19

The cumulative number of daily COVID-19 deaths for each country were downloaded from the JHU CSSE COVID19 Database (date accessed: 19 August 2020) (54). We removed data pertaining to cruise ships, and aggregated data to the country level where data were reported for subnational divisions (e.g. Australia). For countries with overseas territories included in the dataset (e.g. United Kingdom), we excluded the cumulative death counts in those overseas territories. For each country we computed a time series of the daily number of deaths by taking the difference in the cumulative number on consecutive days. When this difference was negative, for example when corrections in the cumulative number were not propagated backwards, we set the value to zero. A relevant outlier in these time series is the addition of 1290 deaths in China on 17 April 2020, while on the days before and after no deaths were recorded. To account for these deaths, we uniformly distributed these deaths over the previous 85 days described by the epidemiological data.

Population data

Country population size estimates were downloaded from the UN Department of Economic and Social Affairs website (https://population.un.org/wpp/Download/Standard/Population/), using the Medium fertility projection for 2020 (55).

Travel and mobility data

To investigate temporal trends in SARS-CoV-2 importation intensity we sought information on the number of travellers entering the UK from each other country for the period from 1 January to 30 April 2020. Incoming travellers comprised both British nationals and resident and visiting citizens of other countries. Estimates were obtained by combining multiple data sources. First, the UK Home Office has provided statistics that describe the number of inbound travellers arriving in the UK by air on each day during this period (https://assets.publishing.service.gov.uk/government/uploads/system/uploads/attachment_data/file/887655/statisticsrelating-to-covid-19-and-the-immigration-system-tables-may-2020-arrivals.ods). This data set provides the daily number of incoming air passengers but not their source country. Second, we obtained the number of tickets sold for inbound flight journeys to the UK along with their origin location from the IATA (for passengers that transfer, the source location is the country from where the whole journey started). We used these numbers to calculate the percentage of arrivals from each country on a monthly basis from January to April 2020. We multiplied the monthly distribution of source destination by the total number of air passenger arrivals in the UK each day to estimate the number of arrivals from each country. Third, we augmented the above air passenger numbers with estimated numbers of incoming travellers arriving per day by short-sea ferry and through the Channel Tunnel (French: Le tunnel sous la Manche). Numbers of short-sea ferry passengers from France, Netherlands and the Republic of Ireland were estimated from monthly statistics obtained from the UK Department of Transport (https://assets.publishing.service.gov.uk/government/uploads/system/uploads/attachment_data/file/908445/spas0107.ods). Within that data set, values are provided for the Republic of Ireland and for “Other EU countries”. The latter total was broken down by country using data from 2019 showing that 72.7% of UK short-sea journeys are with France, 13.6% with the Republic of Ireland, 10.4% with the Netherlands, and 3.3% with other countries (https://www.gov.uk/government/statistics/sea-passenger-statistics-2019-short-sea-routes). Eurotunnel Shuttle vehicle movements from France were obtained from publicly available monthly records (https://www.eurotunnelfreight.com/uk/2020/02/shuttle-traffic-for-january-2020). In the absence of other information we assumed (i) inbound and outbound vehicle movements via the Eurotunnel Shuttle services were equally frequent and (ii) one passenger per truck and 1.5 passengers per passenger vehicle. Inbound Eurostar rail passenger numbers from France and Belgium were estimated from available data and adjusted as far as possible for post-pandemic reduction in travel. Specifically, ∼2m passengers travelled by Eurostar in the first quarter of 2020 (https://www.breakingtravelnews.com/news/article/eurostar-passenger-count-slips-by-a-fifth-in-early-2020).

Monthly Eurostar passenger numbers were then calculated by assuming (i) inbound and outbound journeys were equally frequent, (ii) two thirds of inbound Eurostar journeys originated in France and one third in Belgium, in approximate proportion to the ratio of services, and (iii) the proportional decrease in Eurostar travel volumes during March and April 2020 was equal to that observed for vehicle movements via the Eurotunnel Shuttle. Our estimates do not incorporate estimates of movements across the land border between the UK and the Republic of Ireland. This is unlikely to be problematic as the numbers of infections in the Republic of Ireland was relatively low compared to other potential source countries during the time period of interest.

Epidemiological model

We sought an estimate of the number of individuals in each source country who are (i) infected with SARS-CoV-2 and (ii) able to travel to the UK and initiate a transmission chain. In what follows we refer to these individuals as the “potential initiators of a transmission lineage” (PITL). We conservatively assumed that symptomatic individuals cannot initiate a transmission chain in the UK, either through being prevented from travelling or perfect isolation on arrival. Thus, our estimates of daily SARS-CoV-2 prevalence includes only pre-symptomatic and asymptomatic individuals. Asymptomatic individuals are counted among the PITL as those capable of initiating a transmission lineage at any time while they are still infectious. Figure S16 illustrates the ways in which individuals are counted towards the daily PITL and their potential disease outcomes.

We estimated the daily number of PITL by back-extrapolating the time series of daily numbers of deaths due to COVID-19 in each source country. COVID-19 deaths were used instead of confirmed cases, as we are primarily interested in temporal dynamics rather than absolute values, and death counts are believed to be less sensitive to changes in case definition, reporting delays and differences in the level of surveillance among countries and regions. Estimates of the latent period (infection to becoming infectious), incubation period (infection to onset of symptoms), the infectious duration, and the time between symptom onset and death (in fatal cases) were used to estimate the number of infected individuals who would go on to die from COVID-19, in each stage of the disease, on each day (Fig. S16). We then estimated the total number of infected individuals on each day by multiplying with the reciprocal of the infection fatality rate (IFR).

Estimates of the periods defined above were taken from peer-reviewed sources. Specifically, we assumed that the time from acquiring an infection to becoming infectious is 3 days (56) and the time to symptom onset 5 days (2 days after becoming infectious) (57). The infectious period for patients who recover from the disease was assumed to end 5 days after symptom onset (56) while those who die from the disease are assumed to do so 18 days after the onset of symptoms (58). Given the large numbers of deaths we expect that variation in these timings among individuals will be averaged out and is not considered. We further assumed an asymptomatic proportion of 31% (59) and an IFR of 1%, which is broadly consistent with those found in the literature for China, France, and passengers aboard the Diamond Princess (58, 60, 61). These values correspond to our study period, the spring epidemic of COVID-19; more recent estimates of IFR may vary due to changing treatment regimes and other factors. To examine the sensitivity of our results to the asymptomatic proportion we re-ran our analysis with proportions of 0.18 and 0.78 (the range of published estimates; 62, 63), and found that our results were robust over this range (data not shown). As our main results are fully determined by temporal trends in EII and not absolute numbers, they are invariant to the value of the IFR and we did not perform a sensitivity analysis on it. We did not account for changing levels of infectivity among individuals over the course of their infection.

Using the time series of deaths extracted from the JHU CSSE COVID-19 Database (54), as described above, we obtained estimates of the daily number of PITL in 183 countries from 31 December 2019 to 26 July 2020.

Estimated importation intensity

The daily “estimated importation intensity” (EII) of a country is defined as the product of the proportion of individuals in that country who make up the PITL (as described above) on each day, and the number of individuals who travelled from that country to the UK on that day. The former is estimated by dividing our estimate of the total number of individuals who could potentially initiate a lineage (for each day) by the total population of the country (see the Epidemiological model section). The latter corresponds to the total number of arrivals by air, ferry, and rail on that day (see the Travel and mobility data section). To assist in the subsequent use of the EII, we aggregated all countries with low PITL estimates into a single “other” category. The aggregated countries are those that comprised less than 1% of the cumulative total number of cases as of 1 May 2020 (excluding the UK). This left 53 primary source locations. Maximum EII (Fig. S17) was highest for Spain, (which experienced a large, early epidemic that peaked before inbound passenger numbers declined), followed by France (whose later epidemic peak coincided with high but declining international travel).

Importation lag model

We modelled the TMRCA of an observed transmission lineage (the data observation) as the arrival date of the index patient (of that transmission cluster) in the UK, G, plus a lag time, L, until the first transmission event in the lineage revealed by the data. Given the probability that an importation occurs on day g, f_G(g), and the probability of a lag time of j days, f_L(j), the probability of a TMRCA occurring on day k is v_k, is defined by with . TMRCAs and importation dates are assumed to be independent, so the likelihood for all transmission lineages is the product of the corresponding v_k for each lineage.

This model does not account for incomplete sampling of patients from UK transmission lineages. It is likely that the TMRCA of a small transmission lineage is more recent than the first transmission event after the importation and this issue is potentially further exacerbated by non-random sampling of genome sequences from patients in the lineage (64). We therefore expect shorter lag times for bigger transmission lineages. To account for this sizedependence, we model the average importation lag as a function of lineage size. The functional form of this is given by the equation α + β / n, where α corresponds to the minimal average lag time expected under complete sampling of the lineage and β accounts for the increase in lag time as a smaller proportion of sequences are included in the lineage.

We applied this model to the TMRCA estimates of individual transmission lineages and their sizes as obtained from the MCC trees (see the Identifying transmission lineages section). Values for α and β were found by numerically optimising the likelihood function using random draws from an exponential distribution as initial parameter values. The optimisation procedure was repeated several times to ensure that the algorithm did not become stuck in a local optimum. We further tested whether lineage size affects the importation lag through a likelihood ratio test (LRT) comparing the above model to a nested model without size dependence (β = 0) and found that the size-dependent model is preferred (, p < 0.001). The maximum likelihood estimates for α and β are 0.72 and 28.91 (Fig. S18), respectively.

Travel advice in the UK

The travel advice issued by the Foreign and Commonwealth Office (FCO) of the United Kingdom pertaining to countries and regions affected by COVID-19 was primarily made available through their website (FCO Travel advice: coronavirus (COVID-19) at https://www.gov.uk/guidance/travel-advice-novel-coronavirus). The number of COVID-19 cases in the UK was available via the government website (Coronavirus cases in the UK at https://www.gov.uk/guidance/coronavirus-covid-19-information-for-the-public). Travel advice was also echoed by various news outlets and other information platforms, such as the Public Health Scotland/NHS Scotland Fit for Travel website (https://www.fitfortravel.nhs.uk/). We collected this information by mining archived FCO sites, manually retrieving HTML files corresponding to updates to the URLs provided above and available at the Internet Archive (https://archive.org/). Files were obtained and examined for all dates when changes to the URL were published (18 updates were published in total between 4 February and 23 May 2020). Furthermore, we compared this advice with the Fit for Travel online resource, collected through a similar approach. Where information was insufficient or unclear, we complemented it with data from news outlets to clarify travel advice, which was the case before February 4, when there was no official travel advice (only notifications for novel coronavirus). We collated all the travel advice information into a single standardised table containing types of advice, dates of implementation and countries or geographic regions covered by the advice. The types of advice included both suggestions against specific types of travel versus all but non-essential travel and the recommended period of self-isolation upon return from specific destinations. All of the changes in travel advice were between February 6 and March 23, when specific self-isolation recommendations applied to the general population and not just returning travellers. A summary of the main changes in the UK travel advice across time (in particular, dates when advice for new countries were issued) is presented in Table S5.