Abstract
New Zealand, a geographically remote Pacific island with easily sealable borders, implemented a nation-wide lockdown of all non-essential services to curb the spread of COVID-19. New Zealand has now effectively eliminated the virus, with low numbers of new cases limited to new arrivals in managed quarantine facilities at the border. Here, we generated 649 SARS-CoV-2 genome sequences from infected patients in New Zealand with samples collected between 26 February and 22 May 2020, representing 56% of all confirmed cases in this time period. Despite its remoteness, the viruses imported into New Zealand represented nearly all of the genomic diversity sequenced from the global virus population. The proportion of D614G variants in the virus spike protein increased over time due to an increase in their importation frequency, rather than selection within New Zealand. These data also helped to quantify the effectiveness of public health interventions. For example, the effective reproductive number, Re, of New Zealand’s largest cluster decreased from 7 to 0.2 within the first week of lockdown. Similarly, only 19% of virus introductions into New Zealand resulted in a transmission lineage of more than one additional case. Most of the cases that resulted in a transmission lineage originated from North America, rather than from Asia where the virus first emerged or from the nearest geographical neighbour, Australia. Genomic data also helped link more infections to a major transmission cluster than through epidemiological data alone, providing probable sources of infections for cases in which the source was unclear. Overall, these results demonstrate the utility of genomic pathogen surveillance to inform public health and disease mitigation.
Main Text
New Zealand is one of a handful of countries that aimed to eliminate coronavirus disease 19 (COVID-19). The disease was declared a global pandemic by the World Health Organisation (WHO) on 11 March 2020. The causative virus, severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2)1, was first identified and reported in China in late December 2019, and is the seventh coronavirus known to infect humans, likely arising through zoonotic transmission from wildlife2. Because of its relatively high case fatality rate3-5, and virus transmission from asymptomatic or pre-symptomatic individuals6–7, SARS-CoV-2 presents a significant public health challenge. Due to its high rate of transmission, morbidity and mortality, SARS-CoV-2 has resulted in world-wide lockdowns, economic collapses and led to healthcare systems being overrun.
Since the publication of the first SARS-CoV-2 genome on 10 January 20208, there has been a substantial global effort to contribute and share genomic data to inform local and international communities about key aspects of the pandemic9. Analyses of genomic data have played an important role in tracking the epidemiology and evolution of the virus, often doing so in real time10, and leading to a greater understanding of COVID-19 outbreaks globally11-15.
New Zealand reported its first case on 26 February 2020 and within a month implemented a stringent, country-wide lockdown of all non-essential services. To investigate the origins, time-scale and duration of virus introductions into New Zealand, the extent and pattern of viral spread across the country, and to quantify the effectiveness of intervention measures, we generated whole genome sequences from 56% of all documented SARS-CoV-2 cases from New Zealand and combined these with detailed epidemiological data.
Between 26 of February and 1 July 2020 there were a total of 1,178 laboratory-confirmed cases and a further 350 probable cases of SARS-CoV-2 in New Zealand (a probable case is defined as a person who has returned a negative laboratory result or could not be tested, but the medical officer of health has assigned the case classification based on exposure history and clinical symptoms). Of these combined laboratory-confirmed and probable cases, 55% were female and 45% were male, with the highest proportion of cases in the 20-29 age group (Table 1). Many cases were linked to overseas travel (37%). Geographic locations in New Zealand with the highest number of reported cases did not necessarily reflect the human population size or density in that region, with the highest incidence reported in the Southern District Health Board (DHB) region rather than in highly populated cities (Figure 1). The number of laboratory-confirmed cases peaked on 26 March 2020, the day after New Zealand instigated an Alert Level 4 lockdown – the most stringent level, ceasing all non-essential services and stipulating that the entire population self-isolate (Figure 1). From 23 May 2020, New Zealand experienced 25 consecutive days with no new reported cases until 16 June, when new infections, linked to overseas travel, were diagnosed. All subsequent new cases have been from patients in managed quarantine facilities.
Demographic data for confirmed (n=1178) and probable (n=350) cases of SARS-CoV-2 in New Zealand between 26 February and 1 July 2020. The percentage of genomes sequenced in each category is shown.
(a) Number of laboratory-confirmed cases by reported date, both locally acquired (grey) and linked to overseas travel (blue) in New Zealand, highlighting the timing of public health alert levels 1-4 (‘eliminate’, ‘restrict’, ‘reduce’, ‘prepare’) and national border closures. The number of genomes sequenced in this study is shown over time. (b) Map of New Zealand’s District Health Boards shaded by the incidence of laboratory-confirmed cases of COVID-19 per 100,000 people. (c) Number of laboratory-confirmed cases per District Health Board (DHB) versus number of genomes sequenced, indicating Spearman’s ρ, where asterisks indicate statistical significance (p<0.0001).
We sequenced a total of 649 virus genomes from samples taken between 26 February (first reported case) and 22 May 2020 (to date, the last confirmed case unassociated with managed quarantine facilities). This represented 56% of all New Zealand’s confirmed cases sampled during that period. The data generated originated from the 20 DHBs from across New Zealand. DHBs submitted between 0.1% and 81 % of their positive samples to the Institute of Environmental Science and Research (ESR), Wellington, for sequencing. Despite this disparity, a strong nationwide spatial representation was achieved (Figure 1).
Notably, the genomic diversity of SARS-CoV-2 sequences sampled in New Zealand represented nearly all of the genomic diversity present in the global viral population, with nine second-level A and B lineages from a recently proposed global SARS-CoV-2 genomic nomenclature16 identified. This high degree of genomic diversity was observed throughout the country (Figure 2). The SARS-CoV-2 genomes sampled in New Zealand comprised 24% aspartic acid (SD614) and 73% glycine (SG614) at residue 614 in the spike protein (Figure 2). Preliminary studies suggest that the D614G mutation can enhance viral infectivity in cell culture17. Nevertheless, it is noteworthy that the increase in glycine in New Zealand samples is due to multiple importation events of this variant rather than selection for this mutation within New Zealand. We also inferred a weak yet significant temporal signal in the data, reflecting the low mutation rate of SARS-CoV-2, which is consistent with findings reported elsewhere (Figure 2).
(a) Root-to-tip regression analysis of New Zealand (blue) and global (grey) SARS-CoV-2 sequences, with the determination coefficient, r2 (an asterisk indicates statistical significance; p<0.05). (b) Maximum-likelihood time-scaled phylogenetic analysis of 649 viruses sampled from New Zealand (coloured circles) on a background of 1000 randomly subsampled viruses from the globally available data (grey circles). Viruses sampled from New Zealand are colour-coded according to their genomic lineage16. (c) The number of SARS-CoV-2 genomes sampled in New Zealand within each lineage16. (d) The sampling location and proportion of SARS-CoV-2 genomes sampled from each viral genomic lineage is shown on the map of New Zealand. (e) The frequency of D (blue) and G (red) amino acids at residue 614 on the spike protein over time.
Despite the small size of the New Zealand outbreak, there were 277 separate introductions of the virus out of the 649 cases considered. Of these, we estimated that 24% (95% CI: 23-30) led to only one other secondary case (i.e. singleton) while just 19% (95% CI: 15-20) of these introduced cases led to ongoing transmission, forming a transmission lineage (i.e. onward transmission to more than one individual; Figure 3). The remainder (57%) did not lead to a transmission event. New Zealand transmission lineages most often originated in North America, rather than in Asia where the virus first emerged, likely reflecting the high prevalence of the virus in North America during the sampling period. By examining the time of the most recent common ancestor, or TMRCA, of the samples, we found no evidence that the virus was circulating in New Zealand before the first reported case on 26 February. Finally, we found that detection was more efficient (i.e. fewer cases were missed) later in the epidemic in that the detection lag (the duration of time from the first inferred transmission event to the first detected case) declined with the age of transmission lineages (as measured by the time between the present and the TMRCA; Figure 3).
(a) Frequency of transmission lineage size. (b) The number of samples in each transmission lineage as a function of the date at which the transmission lineage was sampled, coloured by the likely origin of each lineage (inferred from epidemiological data). Importation events that led to only one additional case (singletons) are shown in grey over time. (c) Frequency of TMRCA (the time of the most recent common ancestor) of importation events over time. (d) The difference between the TMRCA and the date as which a transmission lineage was detected (i.e. detection lag) as a function of TMRCA. Spearman’s ρ indicates a significant negative relationship (p<0.01).
The largest clusters in New Zealand were often associated with social gatherings such as weddings, hospitality and conferences18. The largest cluster, which comprised lineage B.1.26, most likely originated in the USA according to epidemiological data, and significant local transmission in New Zealand was probably initiated by a superspreading event at a wedding in Southern DHB (geographically the most southern DHB) prior to lockdown. Examining the rate of transmission of this cluster enables us to quantify the effectiveness of the lockdown. Its effective reproductive number, Re, decreased over time from 7 at the beginning of the outbreak (95% credible interval, CI: 3.7-10.7) to 0.2 (95% CI: 0.1-0.4) by the end of March (Figure 4). The sampling proportion of this cluster, a key parameter of the model, had a mean of 0.75 (95% CI: 0.4-1), suggesting sequencing captured the majority of cases in this outbreak. In addition, analysis of genomic data has linked five additional cases to this cluster that were not identified in the initial epidemiological investigation, highlighting the added value of genomic analysis. This cluster, seeded by a single-superspreading event that resulted in New Zealand’s largest chain of transmission, illustrates the link between micro-scale transmission to nation-wide spread (Figure 4).
Maximum clade credibility phylogenetic tree of New Zealand’s largest cluster with an infection that most likely originated in the USA. Estimates of the effective reproductive number, Re, are shown in violin plots superimposed onto the tree, grouping the New Zealand samples into two time-intervals as determined by the model. Black horizontal lines indicate the mean Re. Tips are coloured by the reporting District Health Board and their location shown on the map.
The dramatic decrease in Re of this large cluster coupled with the relatively low number of virus introductions that resulted in a transmission lineage suggests that implementing a strict and early lockdown in New Zealand rapidly reduced multiple chains of virus transmission. As New Zealand continues to remain free of COVID-19 community transmission, but with positive cases still detected amongst individuals quarantined at the border reflecting high virus incidence in other localities, it is imperative that ongoing genomic surveillance is an integral part of the national response to monitor any re-emergence of the virus, particularly when border restrictions might eventually be eased.
Methods
Ethics statement. Nasopharyngeal samples testing positive for SARS-CoV-2 by real-time polymerase chain reaction (RT-PCR) were obtained from public health medical diagnostics laboratories located throughout New Zealand. All samples were de-identified before receipt by the researchers. Under contract for the Ministry of Health, ESR has approval to conduct genomic sequencing for surveillance of notifiable diseases.
Genomic sequencing of SARS-CoV-2. A total of 733 laboratory-confirmed samples of SARS-CoV-2 were received by ESR for whole genome sequencing. Viral extracts were prepared from respiratory tract samples where SARS-CoV-2 was detected by RT-PCR using WHO recommended primers and probes targeting the E and N gene. Extracted RNA from SARS-CoV-2 positive samples were subject to whole genome sequencing following the ARTIC network protocol (V1 and V3) and the New South Wales (NSW) primer set15.
Briefly, three different tiling amplicon designs were used to amplify viral cDNA prepared with SuperScript IV. Sequence libraries were then constructed using Illumina Nextera XT for the NSW primer set or the Oxford Nanopore ligation sequencing kit for the ARTIC protocol. Libraries were sequenced using Illumina NextSeq chemistry or R9.4.1 MinION flow cells, respectively. Near-complete (>90% recovered) viral genomes were subsequently assembled through reference mapping. Steps included in the pipeline are described in detail online (https://github.com/ESR-NZ/NZ_SARS-CoV-2_genomics).
The reads generated with Nanopore sequencing using ARTIC primer sets (V1 and V3) were mapped and assembled using the ARTIC bioinformatics medaka pipeline (v 1.1.0)19. For the NSW primer set, raw reads were quality and adapter trimmed using trimmomatic (v 0.36)20. Trimmed paired reads were mapped to a reference using the Burrows-Wheeler Alignment tool21. Primer sequences were masked using iVar (v 1.2)22. Duplicated reads were marked using Picard (v 2.10.10)23 and not used for SNP calling or depth calculation. Single nucleotide polymorphisms (SNPs) were called using bcftools mpileup (v 1.9)24. SNPs were quality trimmed using vcflib (v 1.0.0)25 requiring 20x depth and overall quality of 30. Positions that were less than 20x were masked to N in the final consensus genome. Positions with an alternative allele frequency between 20% to 79% were also masked to N. In total, 649 sequences passed our quality control (BioProject: PRJNA648792; a list of genomes and their sequencing methods are provided in Supplementary Table 1).
Phylogenetic analysis of SARS-CoV-2
SARS-CoV-2 sequences from New Zealand, together with 1,000 genomes uniformly sampled at random from the global population from the ~50,000 available sequences from GISAID26 (June 2020), were aligned using MAFFT(v 7)27 using the FFT-NS-2 algorithm. A maximum likelihood phylogenetic tree was estimated using IQ-TREE (v 1.6.8)28, utilising the Hasegawa-Kishino-Yano (HKY+Γ)29 nucleotide substitution model with a gamma distributed rate variation among sites (the best fit model was determined by ModelFinder30), and branch support assessment using the ultrafast bootstrap method31. We regressed root-to-tip genetic divergence against sampling dates to investigate the evolutionary tempo of our SARS-CoV-2 samples using TempEst (v 1.5.3)32. Lineages were assigned according to the proposed nomenclature16 using pangolin (https://github.com/hCoV-2019/pangolin). To depict virus evolution in time, we used Least Squares Dating33 to estimate a time-scaled phylogenetic tree using the day of sampling. With the full set of New Zealand sequences, we used a time-aware coalescent Bayesian exponential growth model available in BEAST (v 1.10.4)34. The HKY+Γ model of nucleotide substitution was again used along with a strict molecular clock. Because the data did not display a strong temporal signal, we used an informative prior reflecting recent estimates for the substitution rate of SARS-CoV-235. The clock rate had a Γ prior distribution as a prior with a mean of 0.8 × 10-3 subs/site/year and standard deviation of 5 × 10-4 (parameterised using the shape and rate of the Γ distribution). Parameters were estimated using Bayesian Markov Chain Monte Carlo (MCMC) framework, with 2 × 108 steps-long chains, sampling every 1 × 105 steps and removing the initial 10% as burn-in. Sufficient sampling was assessed using Tracer (v 1.7.1)36, by verifying that every parameter had effective sampling sizes above 200. Virus sequences were annotated as ‘imported’ (including country of origin) or ‘locally acquired’, according to epidemiological data provided by EpiSurv37. From a set of 1,000 posterior trees, we estimated a number of statistics using NELSI38. We determined the number of introductions of the virus into New Zealand as well as the changing number of local transmission lineages through time, with the latter defined as two or more New Zealand SARS-CoV-2 cases that descend from a shared introduction event of the virus into New Zealand39. Importation events that led to only a single case rather than a transmission lineage are referred to as ‘singletons’. For each transmission lineage and singleton, we inferred the TMRCA.
To estimate Re through time we analysed New Zealand sequences from the clade identified to be associated with a wedding. We used a Bayesian birth-death skyline model using BEAST (v 2.5)40, estimating Re for two time-intervals, as determined by the model, and with the same parameter settings as above. We assumed an infectious period of 10 days, which is consistent with global epidemiological estimates41.
Data Availability
All genomic data is available on GenBank as noted in the manuscript.
Online Supplementary Material
Supplementary Table 1. A list of genomes and which amplification and sequencing method was used in for each case.
Acknowledgements
This work was funded by the Ministry of Health of New Zealand, New Zealand Ministry of Business, Innovation and Employment COVID-19 Innovation Acceleration Fund (CIAF-0470), ESR Strategic Innovation Fund and the New Zealand Health Research Council (20/1018). We thank the ATRIC network for making their protocols and tools openly available and specifically Josh Quick for sending the initial V1 and V3 amplification primers. We thank Genomics Aotearoa for their support. We thank the diagnostic laboratories that performed the initial RT-PCRs and referred samples for sequencing as well as the public health units for providing epidemiological data. We thank the Nextstrain team for their support and timely global and local analysis. We thank all those who have contributed SARS-CoV-2 sequences to GenBank and GISAID databases.
References
