Phylodynamics reveals the role of human travel and contact tracing in controlling COVID-19 in four island nations

Most populated corners of the planet have been exposed to SARS-CoV-2, the coronavirus behind the COVID-19 pandemic. We examined the progression of COVID-19 in four island nations that fared well over the first three months of the pandemic: New Zealand, Australia, Iceland, and Taiwan. Using Bayesian phylodynamic methods, we estimated the effective reproduction number of COVID-19 in the four islands as 1-1.4 during early stages of the pandemic, and show that it declined below 1 as human movement was restricted. Our reconstruction of COVID-19's phylogenetic history indicated that this disease was introduced many times into each island, and that introductions slowed down markedly when the borders closed. Finally, we found that New Zealand clusters identified via standard health surveillance largely agreed with those defined by genomic data. Our findings can assist public health decisions in countries with circulating SARS-CoV-2, and support efforts to mitigate any second waves or future epidemics.


Introduction
Many continued hoping that the epidemic would soon die out and they and their families be spared. Thus they felt under no obligation to make any change in their habits, as yet. Plague was an unwelcome visitant, bound to take its leave one day as unexpectedly as it had come.
Albert Camus ("The Plague", 1947) The respiratory tract illness caused by Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) known as Coronavirus Disease 2019 (COVID-19; Rothan and Byrareddy, 2020) has caused substantial global health and economic impact, and the global epidemic appears far from over. COVID-19 was first detected in the city of Wuhan, Hubei province (China) in late 2019; (Zhu et al., 2020). By 11 March 2020 there were 100,000 cases across 114 countries leading to the World Health Organization (WHO) declaring a pandemic (World Health Organization, 2020b). By the end of April, 3 million cases and 200,000 deaths had been confirmed (World Health Organization, 2020a), with over half of the world population under stay-at-home orders (Sandford, 2020).
COVID-19 has now been characterised extensively. The average number of people each case goes on to infect (known as the effective reproduction number, R e ), for example, has been estimated as 1-6 during the early stages of the outbreak (Liu et al., 2020;Binny et al., 2020;Alimohamadi et al., 2020). Because COVID-19 often causes mild symptoms (Day, 2020b,a), the difference between the true and reported number of cases could be of an order of magnitude (Li et al., 2020a;Lu et al., 2020). Fortunately, R e has since declined in many locations to less than 1 as result of social distancing, travel restrictions, contact tracing, among other mitigation efforts (Binny et al., 2020;Chinazzi et al., 2020).
At the moment of writing, different parts of the world tell contrasting stories about their struggle against COVID-19. New Zealand, Iceland, and Taiwan, for example, have fared comparatively well with smaller infection peaks and lower excess mortalities compared to many mainland countries (World Health Organization, 2020a;FT Visual & Data Journalism team, 2020). Australia was close to elimination, but in late June a second uncontrolled outbreak began in the state of Victoria. Nevertheless the initial experience of these four countries was similar. They all responded to COVID-19 by closing their borders and requiring isolation or quarantine of those who were allowed through (Hale et al., 2020), and either restricting or discouraging human movement (Apple, 2020). These countries provide excellent case studies to characterise the effects of human travel on epidemic spread, a critical endeavor when facing the threat of a COVID-19 second wave and global endemicity.
Mathematical methods for studying epidemics are as diverse as the diseases they target, making varied 2 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted August 6, 2020. . https://doi.org/10.1101 use of surveillance and genetic data, and ranging from models of mathematical epidemiology (Grassly and Fraser, 2008;Binny et al., 2020;Ferretti et al., 2020;Hethcote, 2000) to phylodynamic models (Volz et al., 2013). The latter enjoy the comparable timescales of epidemiological and evolutionary processes, which result from short viral generation times and high viral mutation rates (20-40 mutations per year for SARS-CoV-2; Lai et al., 2020;Li et al., 2020b;Rambaut, 2020). By modelling the transmission process as phylogenetic trees (Fig. 1), phylodynamic and molecular evolution models can be coupled to reconstruct the mutational history of viral genomes (Li et al., 2020b;Rambaut, 2020;Lai et al., 2020), elucidating the epidemiological processes that shaped them.
We have sequenced 217 SARS-CoV-2 genomes from New Zealand, and evaluated the effects of travel and movement restrictions on disease progression in New Zealand, Australia, Iceland, and Taiwan. Finally, we examined the spatiotemporal structure of New Zealand's COVID-19 epidemic, by comparing genetically identified clusters to those independently identified by public health surveillance.

Phylodynamic models
In the study of disease phylodynamics, a complete epidemic can be modelled as a rooted, binary phylogenetic trees whose N tips are all infected individuals (Fig. 1), and whose characteristics and drivers we wish to uncover. In a rapidly growing epidemic the age of this tree is a good proxy for the duration of the infection in the population. Because it is seldom possible to sample the entire infected population, one instead reconstructs a subtree T from a sample of cases of size n, where n < N. This sample provides n viral genomic sequences we refer to as D, as well as each sequence's sampling date, collectively represented by vector y. D and y thus constitute the input data, and are herein referred to as the "samples". In the case of structured phylodynamic models, samples in y are further annotated by "deme" -a geographical unit such as a city, a country, or a continent. A Bayesian phylodynamic model is characterised by the following unnormalised posterior probability density function: where P(D|T , µ c , θ s ) is the phylogenetic likelihood of tree T given alignment D, whose sequences accumulate substitutions at rate µ c under a substitution model with parameters θ s . f (T |θ τ ) is the probability density function of the tree prior sampling distribution, given tree prior parameters θ τ (in this study we 3 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 6, 2020. . https://doi.org/10.1101/2020.08.04.20168518 doi: medRxiv preprint consider a range of phylodynamic models that differ in their tree prior; see below and Supplementary materials). Finally, f pr (µ c , θ s , θ τ ) is the prior distribution on the remaining parameters.

Multitype birth-death model
The multitype birth-death (MTBD) model (Kühnert et al., 2016;Sciré et al., 2020) is a structured phylodynamic model that assumes infected ("I"; Fig. 1) hosts of viral lineages are organised into a set of demes ("types"), and migrate between them at rate m ( Fig. 1). Infected individuals transmit the infection at birth rate λ, and can be sampled ("S"; Fig. 1) and sequenced at sampling rate ψ, after which they are removed ("R"; Fig. 1) with probability r. We set r = 1 by assuming (sampled) confirmed cases are isolated and are unlikely to cause further infections. Removal without sampling happens at death rate µ (Fig. 1).
Here, "death" refers to anything causing hosts to no longer be infectious, including recovery, host death or behavioural changes (e.g., self-isolation).
The probability density function of T under the MTBD sampling distribution (Kühnert et al., 2016) and its parameters θ τ is given by: where π G represents the frequency of lineages in each deme at the root, and O is tree T 's age.

4
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 6, 2020. . https://doi.org/10.1101/2020.08.04.20168518 doi: medRxiv preprint The MTBD model is reparameterised into a computationally and epidemiologically convenient form, with new parameters R e = λ µ+ψr , b = µ + ψr and s = ψ ψ+µ . R e is the disease's effective reproductive number, which is closely related to the basic reproduction number (Stadler et al., 2013), with both being equal at the onset of the epidemic. b is the rate at which individuals become non-infectious, and s is the probability of an individual being a sample (out of all removed individuals). The phylodynamic tree prior under this parameterisation is: Under our MTBD model, R e , b, s, and m are estimated across a series of time intervals (and held constant within intervals) defined by time boundaries t, with t 0 < t 1 < . . . < t f . t is not estimated here (see below andS upplementary material for more details), but fixed at epidemiologically informed dates, or dates derived from a human mobility model; t 0 is the beginning of the infection, and t f is both the last sampling time and the end of the last interval ( Fig. 1). Time intervals are both parameter-and deme-specific.
Finally, we consider four independent geographical models, each characterised by two demes: a target island deme (New Zealand, Australia, Iceland, or Taiwan; referred to as the "I S" deme) and a mutually exclusive "rest-of-the-world" deme ("RW"), which contributes samples from other countries in the world.
All RW parameters are considered nuisance parameters.

Alternative phylodynamic models
Phylodynamic analyses can be carried out with a range of different models. In order to determine how robust our estimates are to model choice, we analysed our SARS-CoV-2 genomic datasets (see below) with three models in addition to MTBD: (i) a simple discrete phylogeography model that treats demes in the same fashion as molecular sites ("DPG", Lemey et al., 2009), (ii) a discrete phylogeography model with two time intervals ("DPG2"), and (iii) a structured coalescent model ("SC", Müller et al., 2018). A brief description of these models and their results can be found in the Supplementary materials.

5
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

SARS-CoV-2 genomic data and subsampling
We sequenced and assembled 217 New Zealand SARS-CoV-2 genomes that are now available on GISAID (Shu and McCauley, 2017), which provided the remaining 3000+ genome sequences analysed here. Due to computational limitations, we subsampled GISAID sequences as follows: (i) two different subsample sizes ("small" vs. "large"), with 200 and 500 RW sequences, and (ii) two sample selection methods ("time" vs. "active") that chose sequences uniformly across time, or proportionally to the number of active cases over time. Up to 250 sequences from the IS deme are included in each alignment (where more than 250 were available, we subsampled them using the respective sampling method).
Permutating the two options above yields four subsampling protocols (detailed in the Supplementary materials), which were applied to the four IS demes, generating a total of 16 SARS-CoV-2 alignments (30 kb long, 0.3% ambiguous sites on average). Each alignment consists of sequences from one of the four islands (I S) and from the rest of the world (RW).

Modelling human movement as tracked by mobile phone data
With the goal of characterising the decrease in human movement over the weeks our SARS-CoV-2 samples were collected, we employed a Bayesian hierarchical sigmoidal model for the analysis of mobile phone data (Apple, 2020;Fig. 2). This analysis revealed that human movement generally decreased in all four islands, with Taiwan being the first to show a declining trend among the islands, but also the least pronounced one. New Zealand, Australia, and Iceland underwent a marked decline in human movement during March 2020, with a slight recovery starting around mid-April.
This analysis informed the three epoch model for R e . The first interval for R e ranged from the origin t 0 to the first reported case in each IS at t 1 . The second describes normal human mobility and ranged from t 1 to the date of mobility reduction t 2 (shown for each IS in Fig. 2). The third ranged from t 2 to the final sample at t f and describes a period of reduced human movement.

SARS-CoV-2 phylodynamics
Inference under the MTBD and SC models suggests they are fairly robust to subsampling protocols, as opposed to DPG and DPG2, which yielded unrealistically small estimates of infection duration (i.e., tree age) and of viral introduction counts from the "active" method. Such model behavior is unsurprising and has been reported elsewhere (De Maio et al., 2015), although more realistic behavior from DPG and DPG2 6 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 6, 2020. was observed when "large" subsamples were used. Analyses under MTBD and SC were prohibitively time-demanding using the larger datasets. For consistency, we only present results from the "small-time" dataset here (but further results can be found in the Supplementary materials).
We estimated R e as 1-2 during the early infection of each respective country and later at less than 1 after human mobility decreased (see credible intervals in Fig. 3). These estimates were not sensitive to sampling strategies. Each island IS saw multiple introductions of SARS-CoV-2 from the rest of the world RW ( Fig. 4 and 5). Under MTBD, we estimated the mean number of sample sequence introductions to be 63, 87, 49, and 65 for New Zealand, Australia, Iceland, and Taiwan, respectively. Conversely, we estimated just one virus export event in each country. Similar results were observed using the SC model, with neither SC nor MTBD producing drastically different estimates between subsampling strategies.
We estimated the substitution rate of SARS-CoV-2 as approximately 9.4 × 10 −4 substitutions per site per year, across the four phylodynamic models, the four islands, and the "small-time" subsampling method.
These 16 (mean) estimates have a range of (6.1 × 10 −4 , 12 × 10 −4 ), and are consistent with other estimates (Lai et al., 2020;Li et al., 2020b;Rambaut, 2020). By following the same pooling procedure, our estimates for the first transmission event range between 30 September and 25 December, 2019. 7 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

Comparison with contact tracing in New Zealand
The New Zealand Ministry of Health (NZMH) defines a significant COVID-19 cluster to comprise ten or more cases connected through transmission, and who are not all part of the same household (New Zealand Ministry of Health, 2020). Five of these clusters are represented by five or more of our samples (Fig. 4a and   c), with a sixth non-significant cluster.
The most stringent definition of a phylogenetic cluster is a group of reciprocally monophyletic samples sharing the same outbreak identifying label. By this standard, only one of the six epidemiological clusters (i.e., six different labels), Matamata, was recovered from genetic data (Fig. 4).
Less stringently, one can define a phylogenetic cluster as a clade that contains all samples sharing a label, with only a single introduction into a focal deme being invoked. Under this definition (i) case exports are allowed, (ii) isolated labelled tips outside of the clade are considered type-1 contract tracing mismatches, and (iii) unlabelled tips that are both from the focal deme and placed within the clade are considered type-2 mismatches. Using this definition, all clusters identified by NZMH were recovered (among those represented by our samples). However, there were five type-1 mismatches: two from an Auckland-based cluster, one from Christchurch, one from Queenstown, and one from Bluff (shown by arrows in Fig. 4). We estimated 18 type-2 mismatches.

Restricting human movement constrains COVID-19 spread
We explored the introduction and transmission history of COVID-19 in New Zealand, Australia, Iceland, and Taiwan -geographically isolated countries that were initially successful in preventing this disease from wreaking havoc as extensively it did in other parts of the world. Our results show that SARS-CoV-2 was introduced into each island on many independent occasions, corroborating the available epidemiological data from health officials (Supplementary materials). Viral imports rapidly declined after the respective islands closed their borders around 20 March 2020 ( Fig. 4 and 5). Closing borders played a significant role in limiting the spread of COVID-19.
Mobile phone data revealed idiosyncratic human movement patterns (Fig. 2). During the early stages of the pandemic, human mobility rapidly declined in both New Zealand and Australia as a result of stay-at-home orders (Hale et al., 2020), as well as in Iceland despite no such official restraints. Taiwan, on the other hand, did not exhibit as much of a decline in human movement, nor did it enforce stay-at-home orders on the general population. A goal of this study was to characterise the extent to which this decline 8 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
We estimated effective reproduction numbers, R e , as approximately 1-2 at the onset of the pandemic in all four islands. These estimates are not as high as in other countries where the pandemic was more widespread and where R e was typically estimated in the range of 2-6 (Liu et al., 2020;Binny et al., 2020;Alimohamadi et al., 2020). Our low estimates reflect the fact that a large fraction of cases in these countries were travel-related imports, and that community transmission rates were low. R e declined as human movement decreased in each of the countries studied, reaching values below 1.0 (Fig. 3).
The case of Taiwan illustrates that an extreme restriction of domestic human movement -officially mandated or otherwise -is not necessarily required for suppressing a pandemic like COVID-19. Instead, extensive contact tracing and the widespread use of face masks in Taiwan were likely key factors . However, a stringent response against COVID-19 greatly restricting mobility such as that adopted by New Zealand is likely to have contributed to the elimination of SARS-CoV-2 by May 9 (Cousins, 2020). New Zealand was the first country with over 1,000 confirmed cases to return to zero active cases, allowing it to re-open its economy, albeit with a still closed border. Most importantly, human morbidity and mortality were limited. Conversely, while Australia also saw a decline in human mobility and initially enjoyed some measure of success against COVID-19, its more lenient approach might underlie this country's ongoing struggle to fully eliminate SARS-CoV-2, as exemplified by the emergence of a second community-transmission led epidemic in the state of Victoria (Dexter, 2020).
Positive outcomes against this and similar pandemics will be dependent on a plethora of factors that must include human movement, but also population density, enforcement of social distancing, hygienic practices such as the use of face masks, the efficiency of contact tracing, high testing rates, and timely border closure. However, none of these should be seen as a silver bullet.

Phylodynamic modelling complements contact tracing
Identifying cases and the sources of their infections is an essential exercise in slowing infectious diseases such as COVID-19 (Hellewell et al., 2020), as recognised by the Director-General of the WHO in his pronouncement: "You cannot fight a fire blindfolded. And we cannot stop this pandemic if we don't know who is infected". Modern health surveillance strategies include medical databases containing travel history , the development of mobile phone applications (Ferretti et al., 2020), and the identification of outbreak subgroups known as "clusters", whose cases share the same origin (New Zealand Ministry of Health, 2020). At the core of health surveillance practice lies contact tracing, the goal of which is to pinpoint the origin of infections by reconstructing the travel history of confirmed cases, and identifying all 9 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 6, 2020. . https://doi.org/10.1101 individuals with whom cases have come in contact.
Nonetheless, contact tracing has its limitations. Establishing the link between cases is restricted, for example, by the ability of cases to recall their recent contacts and travel history, which can lead to type-2 mismatches between epidemiological and genetic clusters. Furthermore, a case may have attended an event causally linked to a cluster but acquired their infection elsewhere and did not infect anyone at the event (type-1 mismatches). The implications of these types of error are manifold. If cluster sizes are underestimated, then so too is the rate of disease spread. If the extent of import-related cases is overestimated, then the impact of international air travel and the extent of community transmission cannot be fully accounted for. Indeed, our results suggest there are an additional 18 unclassified cases in New Zealand that could have been linked to known clusters but were not, and 5 cases that were linked to a cluster where they did not acquire or transmit the infection. However, by and large, our phylodynamic analysis is in agreement with conclusions reached by NZMH.
Overall, we have shown that contact tracing has been accurate among the cases considered, and succeeded to a large degree in identifying individuals belonging to the same infection cluster. We have demonstrated how the rapid real-time availability and assessment of viral genomic data can complement and augment health surveillance strategies. Genome sequencing is increasingly affordable and shared global surveillance resources such as NextStrain and GISAID are becoming more complete and useful.
We advocate incorporating modern methods of genetic epidemiology into standard contact tracing and epidemic surveillance in order to ensure COVID-19 and future epidemics leave as promptly as they come.

10
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

11
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 6, 2020.  12 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted August 6, 2020. . https://doi.org/10.1101/2020.08.04.20168518 doi: medRxiv preprint Clade support: p ≥ 0.9 0.5 ≤ p < 0.9 t1: First reported case in IS t2: Mobility decrease in IS 13 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted August 6, 2020. . https://doi.org/10.1101/2020.08.04.20168518 doi: medRxiv preprint . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted August 6, 2020. . https://doi.org/10.1101/2020.08.04.20168518 doi: medRxiv preprint . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted August 6, 2020. . https://doi.org/10.1101/2020.08.04.20168518 doi: medRxiv preprint . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted August 6, 2020. . https://doi.org/10.1101/2020.08.04.20168518 doi: medRxiv preprint . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted August 6, 2020. . https://doi.org/10.1101/2020.08.04.20168518 doi: medRxiv preprint