Transmission dynamics of SARS-CoV-2 on the Diamond Princess uncovered using viral genome sequence analysis.

An outbreak of coronavirus disease 2019 (COVID-19) caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) occurred aboard the Diamond Princess cruise ship between her January 20 departure and late February 2020. Here, we used phylodynamic analyses to investigate the transmission dynamics of SARS-CoV-2 during the outbreak. Using a Bayesian coalescent-based method, the estimated mean nucleotide substitution rate of 240 SARS-CoV-2 whole-genome sequences was approximately 7.13 × 10-4 substitutions per site per year. Population dynamics and the effective reproductive number (Re) of SARS-CoV-2 infections were estimated using a Bayesian framework. The estimated origin of the outbreak was January 21, 2020. The infection spread substantially before quarantine on February 5. The Re peaked at 6.06 on February 4 and gradually declined to 1.51, suggesting that transmission continued slowly even after quarantine. These findings highlight the high transmissibility of SARS-CoV-2 and the need for effective measures to control outbreaks in confined settings.

Introduction

Coronavirus disease 2019 (COVID-19) is a respiratory infectious disease caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) (Coronaviridae Study Group of the International Committee on Taxonomy of Viruses, 2020, Huang et al., 2020). The initial case was reported in China in December 2019 (Zhou et al., 2020), and the outbreak spread rapidly around the world resulting in an ongoing pandemic (WHO, 2020b). As of September 27, 2020, the total number of confirmed cases had reached 32.7 million, and 991,000 people had died of COVID-19 (WHO, 2020a).

Coronaviruses are a family of enveloped positive-strand RNA viruses that possess a single-stranded, positive-sense RNA genome comprising 26–32 kilobases (kb). Among them, SARS-CoV-2 belongs to the genus Betacoronavirus and has a genome that is approximately 30 kb in length (Lu et al., 2020). Some betacoronaviruses have caused severe disease outbreaks such as severe acute respiratory syndrome (SARS) and the Middle East respiratory syndrome (MERS). The COVID-19 pandemic is the third major outbreak caused by betacoronaviruses in the past two decades. In contrast to SARS and MERS, a substantial proportion of SARS-CoV-2 infections causing COVID-19 have been asymptomatic (He et al., 2020, Mizumoto et al., 2020, Nishiura et al., 2020). Moreover, since asymptomatic and presymptomatic individuals infected with SARS-CoV-2 can transmit the virus (Bai et al., 2020, Ganyani et al., 2020, Rothe et al., 2020), epidemiological data can be unreliable, causing difficulties in understanding transmission dynamics and controlling the spread of infection.

An outbreak of COVID-19 occurred on board the Diamond Princess cruise ship in February 2020. The Diamond Princess sailed from Yokohama on January 20, 2020, arrived at Hong Kong on January 25, and returned to Yokohama on February 3 with 2666 passengers and 1045 crew members. The presumed COVID-19 index case was reportedly a passenger from Hong Kong who boarded the ship in Yokohama on January 20 and developed a mild dry cough on January 23 (Leung et al., 2020). After the passenger returned to Hong Kong and disembarked on January 25, he became feverish on January 29 and was diagnosed with COVID-19 on February 1. After learning that a former passenger had been diagnosed with COVID-19, the Japanese Ministry of Health, Labour, and Welfare quarantined the ship in Yokohama on February 5. Movement of the passengers was restricted, and self-isolation in cabins was implemented. Symptomatic and high-risk individuals were initially tested for SARS-CoV-2 using reverse-transcription polymerase chain reaction (RT-PCR), and then, asymptomatic passengers were tested after February 11 (National institute of infectious diseases, 2020, Nishiura, 2020, Yamahata and Shibata, 2020). Individuals testing positive were transported from the ship to the quarantine location or hospitals (Yamahata and Shibata, 2020). A total of 712 individuals had tested positive among >3600 tests as of August 2020 (Ministry of Health, 2020).

Understanding the patterns of transmission and the dynamics of an infectious agent, among a susceptible population, is essential for establishing effective strategies to control and prevent the spread of infection, particularly when the pathogenic agent is novel (Gardy and Loman, 2018, Ladner et al., 2019). The COVID-19 outbreak on the Diamond Princess cruise ship might provide valuable samples and information about the transmission dynamics of SARS-CoV-2, particularly in a confined setting. In this case, the presumed index case is known, the number of susceptible individuals was constant, timely preventive measures were deployed, and the passengers were extensively tested using RT-PCR (Liu et al., 2020). Although the number of confirmed cases increased soon after returning to port suggesting that the infection had spread among passengers and crew during the voyage, the temporal dynamics of transmission remained elusive. Here, we analyzed the transmission dynamics of SARS-CoV-2 aboard the Diamond Princess using its viral genome sequence in a Bayesian framework. We estimated the effective population size, as well as basic (R0) and effective (Re) reproductive numbers of SARS-CoV-2 during the outbreak and evaluated changes in transmissibility over time.

Materials and methods

SARS-CoV-2 genomic sequence dataset

All publicly available SARS-CoV-2 genome sequences with clinical information as of August 7, 2020 (n = 78,448), were retrieved from the Global Initiative on Sharing All Influenza Data (GISAID) database (Elbe and Buckland-Merrett, 2017). We included sequences with specific collection date known to the day, sampling location, and host information indicating that they had been collected from human hosts. We retained only complete or nearly complete sequences between 29,600 and 31,000 nt that contained <1% ambiguous nucleotides. Sequences were aligned and trimmed according to the reference sequence (NC_045512.2, 29,903 nt) using MAFFT (v.7.453; https://mafft.cbrc.jp/alignment/software/) with the “-addfragments” option. Finally, we excluded sequences with >20 nt gaps, >20 mismatched nucleotides, or concentrated ambiguity calls (>10 in any 50-nt window) as described (Korber et al., 2020) and manually inspected the sequences using AliView (v.1.26; https://ormbunkar.se/aliview/). All retained sequences were randomly resampled up to two per country, year, and month using nextstrain CLI (v.1.16.5) (Hadfield et al., 2018) to maintain the geographic and temporal diversity of the virus. All available isolates from Wuhan sampled in December 2019 were included in the dataset. The first and last 50 nt of the 5′ and 3′ terminals, respectively (1–50 and 29,854–29,903 nt of NC_045512.2), were trimmed as being error-rich. An outlier sequence was removed based on root-to-tip regression analysis using treetime (v.0.7.6; https://github.com/neherlab/treetime), and the resulting global dataset comprised 240 sequences (Supplementary Table 1). Single nucleotide variants were called using snp-sites (v.2.5.1; https://github.com/sanger-pathogens/snp-sites), VEP (v.101.0; https://asia.ensembl.org/info/docs/tools/vep/index.html), and bcftools (v.1.10.2; http://samtools.github.io/bcftools/bcftools.html). Supplementary Fig. 1 shows the characteristics of single nucleotide variants in the reference dataset.

In total, 74 SARS-CoV-2 sequences from infected persons on the Diamond Princess were identified in the GISAID database as of August 7, 2020. Sequences with multiple ambiguous nucleotides (>1%) and long gaps in the coding region were excluded. Finally, 67 sequences were retained and subsequently analyzed (Supplementary Table 2). All sequences from the ship were collected between February 10 and 17, 2020.

Phylogenetic analyses

The general time-reversible model with invariable sites and γ distribution among site rate variation (GTR+G+I) was selected as the best-fit model for the reference dataset using ModelTest-NG (v.0.1.6; https://github.com/ddarriba/modeltest). Phylogenetic trees were constructed using RAxML (v.8.2.9; https://cme.h-its.org/exelixis/software.html) with 1000 replicates using the maximum likelihood (ML) algorithm. Trees were created using FigTree (v.1.4.4; http://evomics.org/resources/software/molecular-evolution-software/figtree/) and were re-rooted using the earliest sample (IPBCAMS-WH-01, GenBank entry: MT019529.1). Bayesian phylogenetic analysis was conducted using a constant coalescent model and the strict molecular clock model implemented in BEAST (v.1.10.4) (Drummond and Rambaut, 2007).

Substitution rate estimation

The substitution rate of the SARS-CoV-2 sequences was calculated in a Bayesian framework using the Markov chain Monte Carlo (MCMC) algorithm implemented in BEAST (v.1.10.4). The substitution rate of the global dataset (n = 240) was estimated under constant, exponential, and Bayesian skyline plot coalescent models with a strict molecular clock model and an uncorrelated relaxed molecular clock model with a lognormal distribution. GTR+G+I was used as a nucleotide substitution model. Each MCMC run contained 50 million states with sampling every 10,000 states. The MCMC convergence and effective sample size (ESS) were analyzed using Tracer (v.1.7.1) (Drummond and Rambaut, 2007), and ESS values >200 were accepted after burn-in. The best-fitting model was selected based on marginal likelihood estimation using path sampling and stepping-stone sampling methods (Baele et al., 2012). Root-to-tip regression of the ML tree was analyzed using Tempest (v. 1.5.3). The root was determined according to the coefficient of determination, R2.

Phylodynamic analyses

The Bayesian skyline plot is a coalescent-based nonparametric method for estimating past population dynamics over time from molecular sequence data (Drummond et al., 2005). The SARS-CoV-2 population dynamics on the Diamond Princess were inferred using a Bayesian skyline plot coalescent model with the strict molecular clock model implemented in BEAST (v.1.10.4). The GTR+G+I substitution model was used, and the prior for the nucleotide substitution rate was set based on the values obtained in the rate estimation. The MCMC run contained 100 million states, with sampling every 10,000 states. We accepted ESS values >200 after a 10% burn-in, and three independent runs were combined using LogCombiner (v. 1.10.4) implemented in the BEAST package. The effective population sizes were plotted using R (v.4.0.2, R Core Team, 2018; https://www.R-project.org/).

Basic reproductive number estimation

The R0 is the predicted number of secondary cases that one case would generate in a completely susceptible population (Delamater et al., 2019). Generally, an outbreak should continue if R0 has a value >1, but not if R0 is <1. The R0 was estimated using the linear equation, R0 = rD + 1, where r is the exponential growth rate and D is the average duration of infectiousness (Pybus et al., 2001, Wallinga and Lipsitch, 2007). We calculated the exponential growth rate of 67 sequences derived from persons aboard the cruise ship using an exponential coalescent model with the strict molecular clock model implemented in BEAST2 (v.2.6.3). Three independent runs were combined. The average duration of infectiousness was calculated as 9.3 days (95% CI, 7.8–10) (He et al., 2020), and the R0 was calculated using R (v.4.0.2) (Supplementary information). Doubling time was calculated as ln(2)/r. Limitations of this method include biases due to the duration of infectiousness and the estimated substitution rate (Pybus et al., 2001).

Effective reproductive number estimation

The birth–death skyline model is based on a forward-in-time model of transmission, death/recovery, and sampling (Stadler et al., 2013). Temporal changes in Re during the outbreak were estimated using the birth–death skyline contemporary model implemented in BEAST2 (v.2.6.3). The GTR + G + I substitution and strict molecular clock models were used with nucleotide substitution rate estimates. For the birth–death model analysis, the epidemiological parameters comprising the rate of becoming noninfectious (δ), sampling proportion (ρ), Re, and the origin of the epidemic were established based on previous estimates (Stadler et al., 2013, Byrne et al., 2020, He et al., 2020, Lai et al., 2020a). The rate of becoming noninfectious was calculated as the reciprocal of the average duration of infectiousness (9.3 days, 95%; CI 7.8–10) (He et al., 2020), which corresponds to 39.4/year (95%CI: 36.5–46.8). A gamma prior with parameters α = 1.7 and β = 23 (mean, 39.1 year−1; 95% CI, 3.59–116]) was taken as a less informative prior to allow the broader range of the period. Sampling proportion (ρ) is the proportion of lineages sampled in the analysis. The number of infected individuals up to February 17 was expected to range between 454 and 3711, which corresponds to the cumulative cases up to February 17, and the total number of individuals aboard the ship, respectively. Therefore, we calculated ρ to be between 0.018 and 0.14, and set a beta prior with parameters α = 2 and β = 15 (mean, 0.118; 95% CI, 0.0155–0.302). The Re of SARS-CoV-2 has been reported to be around 2–3 and 11 in the community-based infection and confined settings, respectively. A lognormal prior with a log-mean of 1 and an SD of 1.25 (mean, 5.94; 95% CI, 0.235–31.5) was used in this study. The origin of the epidemic is considered as the time point at which the index case became infected. We used 1/X prior as the origin of the epidemic prior. Three independent runs were combined as described previously herein, and Re was plotted using R (v.4.0.2) with bdskytools (https://github.com/laduplessis/bdskytools).

Confirmed cases on the Diamond Princess

The numbers of daily confirmed and cumulative cases obtained from the website of the Japanese Ministry of Health, Labour, and Welfare (https://www.mhlw.go.jp/index.html) and previous studies (Mizumoto et al., 2020, Yamahata and Shibata, 2020) were plotted using R (v.4.0.2).

Results

SARS-CoV-2 sequences and phylogenetic analyses

An ML phylogenetic tree was constructed using the SARS-CoV-2 whole-genome sequences from the 67 persons aboard the Diamond Princess and the 240 global reference sequences (Fig. 1 , Supplementary Fig. 2). Although branch bootstrap values were low, all sequences from the cruise ship clustered with some isolates collected in other countries. The Bayesian phylogenetic analysis also showed a similar clustering for samples from the ship (Supplementary Fig. 3), suggesting that the outbreak originated from a single source.

Fig. 1

Maximum likelihood phylogenetic tree based on SARS-CoV-2 whole-genome sequences of 240 global reference strains and 67 isolates sampled from persons on board the Diamond Princess. Branches of sequences from the Diamond Princess and those from Asian countries are colored red and purple, respectively. The tree was rooted using the earliest sample in the dataset (IPBCAMS-WH-01, GenBank entry: MT019529.1). Branch length reflects units of nucleotide substitutions per site according to the provided scale. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Using the reference sequence dataset, the substitution rate was estimated under several settings (Supplementary Table 3). However, convergence was not achieved using the Bayesian skyline plot model. We selected strict clock and exponential coalescent models as providing the best fit based on marginal likelihood estimations. The estimated mean substitution rate was 7.13 × 10−4 substitutions per site per year (95% highest posterior density [HPD] intervals, 6.16–8.16 × 10−4). The estimated mean time to the most recent common ancestor (tMRCA) was 2019.80 (95% HPD, 2019.71–2019.88), which corresponded to October 19, 2019 (95% HPD, September 18–November 18, 2019). Root-to-tip regression analysis of the ML phylogenetic trees against sampling time was conducted using Tempest (v. 1.5.3). The estimated mean evolutionary rate was 7.82 × 10−4 substitutions per site per year (Supplementary Fig. 4), and the estimated tMRCA was 2019.83 (November 30, 2019).

The SARS-CoV-2 population dynamics on the Diamond Princess were inferred using the Bayesian skyline plot coalescent and strict molecular clock models. The Bayesian skyline plot showed that the effective population size began to increase around January 30 and exponentially increased between February 2 and 6, 2020 (Fig. 2 a). Subsequently, the population growth slowed down but continued on a slightly increasing trajectory, followed by a plateau phase from approximately February 10, 2020. The estimated mean tMRCA was 2020.066 (95% HPD, 2020.027–2020.095), corresponding to January 25, 2020 (95% HPD, January 10–February 4).

Fig. 2

Transmission dynamics of SARS-CoV-2 on board the Diamond Princess. (a) Effective population size of SARS-CoV-2 on board the Diamond Princess over time, (b) effective reproductive number (Re) over time, and (c) daily and cumulative numbers of confirmed cases. Effective population size and Re were estimated from 67 whole-genome sequences sampled from persons on board the Diamond Princess. Thick black line in (a) and (b), median estimate; colored area, 95% highest posterior density (HPD) confidence intervals of the estimate. Dashed (a) and dotted (b) lines, time to the most recent common ancestor and estimated time of outbreak origin, respectively.

Basic and effective reproductive numbers over time

To estimate R0 and the doubling time of infections, the exponential growth rate for 67 sequences from the cruise ship was calculated under an exponential coalescent model with a strict clock model. The estimated exponential growth rate was 61.5 per year (95% HPD, 1.66–143.5), corresponding to 0.168 per day (95% HPD, 0.00454–0.392). The estimated R0 was 2.56 (95% confidence interval [CI], 1.04–4.65), and the direct estimation of the mean doubling time was 4.11 days (95% CI, 1.76–152.4; Table 1 ).

Table 1

Estimated epidemiological parameters based on exponential growth rate.

	Mean estimates	95% CI
Basic reproductive number (R0)	2.56	1.04, 4.65
Doubling time (days)	4.11	1.76, 152.4

We estimated epidemiological parameters and temporal changes in Re for the outbreak using the birth–death skyline contemporary model (Table 2 , Fig. 2b). The number of mean infected persons was calculated as 682 (95% HPD, 293–17000) based on the estimated sampling proportion (ρ). The Bayesian birth–death skyline plot showed a median Re of 2.84 (95% HPD, 0.0204–19.6) at the origin of the epidemic, which we estimated had been between January 21 and 22, 2020. The Re remained at 2–3 during the cruise, reached a high of 6.06 (95% HPD, 0.0379–31.8) on February 4, and then gradually declined to 1.51 (0.0320–6.33) on February 17, reflecting the effects of placing the ship in quarantine.

Table 2

Estimated epidemiological parameters using a birth–death model.

	Mean estimates	95% HPD
Origin of the epidemic (y)	0.0722 (January 21, 2020)	0.0360, 0.124 (Jan 2–Feb 3, 2020)
Rate of becoming noninfectious per year	30.4	4.83, 64.1
Sampling proportion (ρ)	0.0982	3.93 × 10−3, 0.229
Tree root (tMRCA, year)	0.0573	0.0334, 0.0858

HPD, highest posterior density; tMRCA, time to the most recent common ancestor.

Fig. 2c shows daily confirmed and cumulative cases over time. Symptomatic and subsequently asymptomatic passengers and crew were tested using RT-PCR when the ship returned to Yokohama on February 3. The first results were obtained on February 5, 2020. The number of confirmed cases increased substantially, and the cumulative number exponentially increased after February 7, which reflected the estimated exponential increase of infections on the cruise ship after February 2.

Discussion

The present study investigated transmission dynamics of SARS-CoV-2 on the Diamond Princess cruise ship based on viral genome data using a Bayesian framework. Previous studies have estimated these transmission dynamics on the ship based on mathematical models using the daily number of confirmed cases or the daily number of symptomatic confirmed cases (Mizumoto and Chowell, 2020, Nishiura, 2020, Rocklöv et al., 2020). However, these epidemiological data could have been biased because of the underestimation of infections owing to the limited diagnostic test capacity in the early phase of the outbreak and the limited sensitivity of RT-PCR (Woloshin et al., 2020). In addition, the origin of the outbreak might be difficult to estimate due to the lack of diagnostic tests on the voyage. In contrast, recent advances in genomics and computing have facilitated the extraction of valuable information from viral genome data, and analytic methods have been developed to infer transmission history and past population dynamics of the pathogen with limited epidemiological data (Pybus et al., 2001, Drummond et al., 2005, Stadler et al., 2013). We aimed to elucidate the dynamics of transmission from the perspective of genomic epidemiology, which would contribute to further understanding of the outbreak.

SARS-CoV-2 can spread rapidly via close contact and in confined spaces (Li et al., 2020, Rocklöv and Sjödin, 2020), such as cruise ships (Moriarty et al., 2020, Sando et al., 2020, Sekizuka et al., 2020b) and airplanes (Hoehl et al., 2020, Pavli et al., 2020, Speake et al., 2020). Our results showed the transmission dynamics of the SARS-CoV-2 outbreak in a high-transmission environment. The estimated origin of the epidemic was January 21, which coincided with the date when the presumed index case boarded the ship (Leung et al., 2020). The estimated R0 was 2.56 (95% HPD, 1.04–4.65), which was in the range of the community-level R0 (2.2–2.7) (Li et al., 2020, Wu et al., 2020) and in lower range of previous estimates for the ship (2.28–14.8) (Rocklöv et al., 2020, Zhang et al., 2020b). The effective population size began to increase around January 30 and exponentially surged from February 2, before the start of quarantine on February 5. These findings indicate that most infections had already occurred before quarantine, which is consistent with previous findings (Mizumoto et al., 2020, National institute of infectious diseases, 2020, Nishiura et al., 2020). When we compare our findings with epidemiological data, the exponential increase in confirmed cases after February 7 appears to follow the estimated substantial growth in infections by approximately 5 days. This time interval is in close agreement with the incubation period of the virus, which is on average 5.2 days with a 95th percentile of the distribution at 12.5 days (Li et al., 2020). Although the number of confirmed cases might be affected by the limited capacity of RT-PCR tests or change of test strategies, this paralleled increase confirms that our estimates reasonably reflect the increase in infections. Previous studies also used the number of confirmed febrile cases by date of symptom onset (Mizumoto and Chowell, 2020, National institute of infectious diseases, 2020, Nishiura, 2020, Tsuboi et al., 2020). The number of confirmed febrile cases surged on February 7 and gradually declined afterwards. The increase in febrile cases might reflect the estimated exponential increase in infections in the beginning of February, although the number of confirmed febrile cases might be biased by the distribution of thermometer among all individuals and the establishment of the Fever Call Center on February 7. The interval between the estimated increase in infections and the observed peak of daily febrile cases appears to be relatively close, considering the average incubation period of 5.2 (95% CI: 4.1–7.0) days (Li et al., 2020). However, the incubation period of COVID-19 varies across individuals within a population, and factors such as age and disease severity affect the length of the incubation period (Dai et al., 2020, Kong, 2020, Lai et al., 2020b). Therefore, the interval between our estimates and the number of febrile cases could be explained by the temporal change of infected subpopulations with different backgrounds.

The modern cruise ships carry large numbers of people in limited spaces for a certain period, which facilitate the transmission of infectious diseases and outbreaks (Mitruka and Wheeler, 2008, Kak, 2015). The Diamond Princess had many passengers and crew from various countries, and her cruise lasted more than 2 weeks. The facilities on the ship provided semi-closed spaces for passengers and crew, which contributed to high transmission of SARS-CoV-2. Large cocktail parties during the cruise played an important role in the dissemination of infections among passengers and crew (Yamagishi et al., 2020). From the perspective of virology, SARS-CoV-2 has highly contagious characteristics, including a higher R0, longer incubation period, shorter interval between symptom onset and infectivity, and higher proportions of asymptomatic persons and patients with mild illness, compared to those with influenza, SARS-CoV, and MERS-CoV (Petersen et al., 2020). These environmental and virologic factors might contribute to the outbreak of SARS-CoV-2 onboard the cruise ship.

The Re peaked at 6.06 on February 4 and then decreased to 1.51. The effective population size continued to slightly increase after quarantine and reached a plateau around February 10. Although quarantine greatly reduced the number of infections, our results suggested that the infection spread slowly even after quarantine, which is consistent with previous epidemiological studies using mathematical models (Mizumoto and Chowell, 2020, Nishiura, 2020) and haplotype network analysis using viral genome sequences (Sekizuka et al., 2020a). Strict implementation of quarantine on the cruise ship was encumbered by difficulties in developing effective traffic lines or zones between infectious and noninfectious passengers and in providing daily support and medical procedures for passengers (Yamahata and Shibata, 2020). The former possibly promoted transmission among passengers after quarantine through shared onboard facilities or cabins (Sekizuka et al., 2020a). The latter necessitated activities of the crew, which put members at risk of infection. A previous study estimated that the epidemic peaked among passengers and crew between February 2–4 and February 8–10, respectively (Nishiura, 2020). The epidemic period among the crew coincided with the increase in the effective population size after quarantine, suggesting that the slight increase in infections reflected the spread of infection among crew members. Further investigations are essential to understand the transmission dynamics and mechanism of SARS-CoV-2, including the possibility and effect of airborne transmission (Centers for Disease Control and Prevention, 2020). It is also necessary to develop effective strategies for preventing and controlling outbreaks in confined settings. The development of vaccines and anticoronaviral drugs are urgent for preventing future COVID-19 epidemic (Basu et al., 2020, Shin et al., 2020, Zhang et al., 2020a). Based on lessons learnt from the outbreak aboard the Diamond Princess, appropriate protective measures, and separation of clean and contaminated areas should be implemented comprehensively under the supervision of infectious disease specialists. In addition, activities of staff should be restricted to a minimum under quarantine, and suspected infected staff should be isolated as early as possible (Tokuda et al., 2020). Finally, cooperation of various specialists and an established chain of command are essential for controlling a large-scale outbreak in confined settings (Tokuda et al., 2020). Comprehensive and practical protocols should be designed for future epidemic events.

Our study had certain advantages when compared with previous studies on the transmission dynamics of SARS-CoV-2 aboard the Diamond Princess using mathematical models and the number of laboratory-confirmed febrile cases (Mizumoto and Chowell, 2020, Nishiura, 2020). First, our study takes asymptomatic cases into consideration. Although previous studies focused on symptomatic population, a substantial proportion of SARS-CoV-2 infections is asymptomatic and has a significantly longer duration of viral shedding than that in the symptomatic population (Long et al., 2020, Mizumoto et al., 2020, Yanes-Lane et al., 2020), suggesting a critical role of asymptomatic individuals in the outbreak. In contrast, most sequences in this study were obtained from both symptomatic and asymptomatic patients through extensive tests, enabling our study to reflect the entire dynamics of transmission on the ship. Second, our analysis used objective data such as viral genome sequences and sampling dates. In contrast, the dates of fever onset can be affected by factors such as the background of individuals and testing situation on the ship. For instance, older adults with SARS-CoV-2 infection become more symptomatic and have a longer incubation period (Dai et al., 2020, Yamagishi et al., 2020). In addition, distribution of thermometers among all individuals and the establishment of the Fever Call Center on February 7 could increase the detection of symptoms (Tsuboi et al., 2020). The data used in this study are more robust against these potential biases. Third, our study investigates the origin of the epidemic using viral genomic data. Our results showed that the outbreak had a single origin and that the epidemic started immediately after departure. However, in our study, one of the disadvantages was that the association between transmission dynamics and epidemiological factors could not be analyzed. Factors such as age, sex, location of cabin, and membership might largely affect transmission patterns. It was technically difficult to combine this information with phylogenetic analysis in this study. Therefore, we could not analyze the transmission dynamics in each subpopulation or between subpopulations. Second disadvantage is the potential sampling bias and sequencing errors due to next-generation sequencing techniques. It is difficult to obtain a complete and high-quality viral sequence from a sample with a low viral load. The viral load can be affected by the sampling procedure, immunological characteristics of individuals, and interval between infection and sampling. Therefore, there is a possibility that specific groups were excluded from our study.

This study has several potential limitations. First, we used SARS-CoV-2 sequences sampled over a period of 8 months. The short time range and relatively low diversity of SARS-CoV-2 might have affected the estimation of the substitution rate. To mitigate this challenge, we selected sequences to represent the temporal and geographical diversity of SARS-CoV-2 at the time while minimizing the impact of sampling bias. Epidemiological parameters of SARS-CoV-2 infection are currently limited. Although the prior used herein was set based on available epidemiological data, these settings might need modification based on future data. However, our results were consistent with those of previous epidemiological studies using mathematical models. Regarding the application of the Bayesian Skyline plot for epidemic data, sequences should be sampled sparsely from a large population under the assumption of a coalescent model. Sequences used in this study were less than 10% of the total confirmed cases, which might be sufficiently sparse (Fu, 2006). In addition, the HPD interval under a coalescent model tends to be overconfident when the number of lineages in the tree gets close to population size (Stadler et al., 2013). Therefore, the HPD interval in the early phase of an epidemic should be interpreted with caution. Re was estimated using the birth–death skyline contemporary model in this study. The birth–death model can be biased largely when sampling process is mis-specified in the use of serially-sampled sequences (Volz and Frost, 2014). It was difficult to establish a prior reflecting the complex sampling process on the ship, and convergence was not achieved. Because sequences were obtained almost contemporaneously, we used the birth–death contemporary model with the sampling proportion prior that can be set with reasonable certainty. This model selection could affect the estimated change of Re within approximately 1 day; however, our results would be robust against the irregular sampling process on the ship. Finally, we analyzed sequences as they were sampled from only one group. However, the transmission patterns differed between the passengers and crew, which might have affected our findings. The discrepancies in transmission patterns between passengers and crew can be explained by different living patterns on the cruise ship, locations of cabins in distinct decks, and extent of activity after quarantine (Nishiura, 2020, Tsuboi et al., 2020, Yamagishi et al., 2020). The incorporation of the different transmission patterns into phylogenetic analyses would contribute to more accurate projections of transmission dynamics. However, it is difficult to incorporate phylogenetic analysis with complicated epidemiological models, including age, structure of host contact rates, or multiple host populations with different life history traits (Volz et al., 2013). Alternatively, we could analyze sequences from a specific subpopulation, for instance passengers or crew. We identified four sequences from the crew in our dataset, which were too few to analyze transmission dynamics among the passengers and crew. Future studies involving more genetic and finer epidemiological data could facilitate further understanding of the complexity of SARS-CoV-2 transmission in confined spaces.

In conclusion, we investigated the transmission dynamics associated with the COVID-19 outbreak aboard the Diamond Princess cruise ship. Based on the estimated substitution rate, phylodynamic analysis showed that the outbreak on the cruise ship occurred between late January and early February, before quarantine started on February 5, 2020. Although quarantine was clearly effective, the infection continued to slowly spread, even after the countermeasures were initiated. More effective measures should be considered to control outbreaks, particularly in confined settings.

Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

CRediT authorship contribution statement

Kunikazu Hoshino: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Visualization, Writing - original draft. Tatsuji Maeshiro: Writing - review & editing. Nao Nishida: Writing - review & editing. Masaya Sugiyama: Writing - review & editing. Jiro Fujita: Writing - review & editing. Takashi Gojobori: Conceptualization, Investigation, Supervision, Writing - review & editing. Masashi Mizokami: Conceptualization, Investigation, Supervision, Project administration, Writing - review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.