Viral transmission and evolution dynamics of SARS-CoV-2 in shipboard quarantine.

OBJECTIVE: To examine transmission and evolution of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) in shipboard quarantine of the Diamond Princess cruise ship.
METHODS: We obtained the full SARS-CoV-2 genome sequences of 28 samples from the Global Initiative on Sharing All Influenza Data database. The samples were collected between 10 and 25 February 2020 and came for individuals who had been tested for SARS-CoV-2 during the quarantine on the cruise ship. These samples were later sequenced in either Japan or the United States of America. We analysed evolution dynamics of SARS-CoV-2 using computational tools of phylogenetics, natural selection pressure and genetic linkage.
FINDINGS: The SARS-CoV-2 outbreak in the cruise most likely originated from either a single person infected with a virus variant identical to the WIV04 isolates, or simultaneously with another primary case infected with a virus containing the 11083G > T mutation. We identified a total of 24 new viral mutations across 64.2% (18/28) of samples, and the virus evolved into at least five subgroups. Increased positive selection of SARS-CoV-2 were statistically significant during the quarantine (Tajima's D: -2.03, P < 0.01; Fu and Li's D: -2.66, P < 0.01; and Zeng's E: -2.37, P < 0.01). Linkage disequilibrium analysis confirmed that ribonucleic acid (RNA) recombination with the11083G > T mutation also contributed to the increase of mutations among the viral progeny.
CONCLUSION: The findings indicate that the 11083G > T mutation of SARS-CoV-2 spread during shipboard quarantine and arose through de novo RNA recombination under positive selection pressure. (c) 2021 The authors; licensee World Health Organization.

Introduction

On 31 December 2019, Chinese authorities alerted the World Health Organization (WHO) of an outbreak of a novel coronavirus causing severe illness with pneumonia-like symptoms. The virus was later named severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). WHO declared the outbreak a Public Health Emergency of International Concern on 30 January 2020 and as of 17 December 2020, more than 76 million confirmed cases had been reported and 1.6 million people had died.

SARS-CoV-2 contains a single positive stranded RNA (ribonucleic acid) of 30 kilobases, which encodes for 10 genes, and like other RNA viruses, evolves through random nucleotide substitutions. Early sequencing of samples from coronavirus disease 2019 (COVID-19) patients on 30 January 2020 showed only six mutations compared with the first isolate WIV04, and hence suggested that the virus had a short infection history in humans. By April 2020, researchers had found more than 131 mutations in SARS-CoV-2 across the 103 sequenced viral genomes and they estimated that the virus accumulates about one to two mutations per month.

On 25 January 2020, a passenger disembarked the Diamond Princess cruise ship in Hong Kong Special Administrative Region, China, and on 1 February, the passenger tested positive for SARS-CoV-2. When the ship docked in Yokohama, Japan, on February 3, the 3711 cruise passengers and crew members had to quarantine on the ship. As of 21 February 2020, 712 (19.2%) of these individuals had tested positive for SARS-CoV-2, of which 331 (46.5%) were asymptomatic at the time of testing. Among the 381 symptomatic individuals, 37 (9.7%) required intensive care and nine (2.4%) died.,

The shipboard quarantine provided a closed environment to observe the SARS-CoV-2 transmission and adaptation independently from other infectious resources., This environment presents an ideal static population, with little interfering noise, to measure the viral phylodynamics from the COVID-19 outbreak. We therefore decided to use this opportunity to study de novo evolution of SARS-CoV-2 in a closed population.

Methods

Data resources

Viral sequences and sequencing methods are available in the Global Initiative on Sharing All Influenza Data (GISAID) and GenBank® databases. From these databases, we downloaded sequences and annotations of the isolates from the cruise ship, as well as the reference genomes of SARS-CoV-2 isolates PBCAMS-WH-04 (accession no. MT019532), WIV04 (MN996528), Hu-1 (NC045512) and WHU01 (MN98868).

Statistical and phylogenetics analyses

We aligned FASTA files of viral sequences using MAFFT 7 software (Kazutaka Katoh, Research Institute for Microbial Diseases, Osaka, Japan)., To analyse phylogenetic relationships between viral sequences, we used the neighbour-joining method and Jukes–Cantor substitution model with setting bootstrap resampling number as five. We generated the rectangular phylogenetic tree using Archaeopteryx with Java plug-in of MAFFT 7. For radial phylogenetic tree, we first exported the tree file as Newick format by MAFFT; the FigTree software (version 1.4.2, Andrew Rambaut, University of Edinburgh, Edinburgh, United Kingdom of Great Britain and Northern Ireland) was used to transform and display the cladogram. To illustrate the viral subgroups with 11803G > T mutation, we rooted an unrooted tree by introducing the bat SARS-like CoV WIV16 (accession no. KT444582) as an outgroup virus.

Selection pressure

To determine whether the viral genome undergoes neutral or non-neutral evolution, that is, genetic variations of viral genomes are due to randomly genetic drift or under natural selection pressure, we used MEGA7 software to calculate Tajima’s D test of neutrality. This method compares the number of mutations per site with the nucleotide diversity (the mean pairwise difference between sequences). To compare the number of derived singleton site mutations – that is, single base mutations that occur only once in a given population – and the mean pairwise difference between sequences, we calculated Fu and Li’s D, using online PopSc calculator (Shi-Yi Chen, Sichuan Agricultural University, Chengdu, China). We used the same calculator for calculating Zeng’s E, which measures changes in high-frequency variants.– To calculate the P values for Tajima’s D, Fu and Li’s D and Zeng’s E values, we used DnaSP 6 software (Julio Rozas, University of Barcelona, Barcelona, Spain).

To investigate the linkage disequilibrium, that is the non-random assortment of alleles at different loci, of SARS-CoV-2 genomes, we first converted 148 SARS-CoV-2 genomic sequences using SNP_tools plug-in in Excel (Microsoft, Redmond, United States of America, USA) to create a baseline. We downloaded these sequences from GISAID. Using HaploView software, version 4.1 (Broad Institute, Cambridge, USA), we measured and plotted the normalized values (D’) of the coefficient of linkage disequilibrium (D). We obtained D’ by dividing D with D, where D is the theoretical maximum difference between the observed and expected haplotype frequencies. We also calculated the log of the odds of there being a disequilibrium between two loci and the squared coefficient of correlation (r) using the same software. In the absence of evolutionary forces or natural selection, the D’ converges to zero along the time axis at a rate depending on the magnitude of the recombination rate between the two loci. We used the χ test to examine if the obtained linkage disequilibrium was statistically significant. To detect positive RNA recombination, we plotted 95% confidence bounds for D’ using HaploView. Pairs are thought to be in strong linkage disequilibrium if the upper 95% confidence bound is above 0.98 (that is, consistent with no recombination) and the lower bound is above 0.7. Conversely, strong evidence for recombination is defined if pairs for which the upper confidence bound of D’ is less than 0.9. We searched a solid spine of strong linkage disequilibrium running from one marker to another along the legs of the triangle in the linkage disequilibrium chart to determine the haplotype block.

Results

Viral variants

A total of 28 specimens with viral sequences were available for this analysis, including 25 samples from the United States and three samples from Japan. Table 1 lists the characteristics of the viral sequences. Genetic variations of viral sequences were present in 71.4% (20/28) of samples. The sequences of eight samples were completely identical to the Wuhan isolates PBCAMS-WH-04, WIV04, Hu-1 and WHU01. A total of 24 new substitution mutations were identified in 18 samples (Table 1 and Fig. 1).

Table 1

Characteristics of SARS-CoV-2 genomes from samples taken from people on shipboard quarantine, February 10 to February 25, 2020

GISAID name	Accession no.	Specimen source	Collection date	Genomic changea	Type of mutation	Gene/protein	Amino acid change
hCoV-19/USA/CruiseA-1/2020	EPI_ISL_413606	Nasopharyngeal swab	17 Feb	3099C > T	Missense	Orf1ab/NSP3	T945I
				28378G > T	Synonymous	N	NA
hCoV-19/USA/CruiseA-2/2020	EPI_ISL_413607	Nasopharyngeal swab	18 Feb	28409C > T	Missense	N	P46S
				29736C > T	Non-coding	3′-UTR	NA
hCoV-19/USA/CruiseA-3/2020	EPI_ISL_413608	Nasopharyngeal swab	18 Feb	No change	NA	NA	NA
hCoV-19/USA/CruiseA-4/2020	EPI_ISL_413609	Nasopharyngeal swab	21 Feb	1385C > T	Synonymous	Orf1ab/NSP2	NA
				29230C > T	Synonymous	N	NA
				29635C > T	Non-coding	3′-UTR	NA
hCoV-19/USA/CruiseA-5/2020	EPI_ISL_413610	Oropharyngeal swab	21 Feb	No change	NA	NA	NA
hCoV-19/USA/CruiseA-6/2020	EPI_ISL_413611	Nasopharyngeal swab	21 Feb	11410G > A	Synonymous	Orf1ab/NSP6	NA
				26326C > T	Synonymous	E	NA
hCoV-19/USA/CruiseA-7/2020	EPI_ISL_413612	Nasopharyngeal swab	17 Feb	3738C > T	Missense	Orf1ab/NSP3	P1158S
				11410G > A	Synonymous	Orf1ab/NSP6	NA
				26326C > T	Synonymous	E	NA
hCoV-19/USA/CruiseA-8/2020	EPI_ISL_413613	Nasopharyngeal swab	17 Feb	9474C > T	Missense	Orf1ab/NSP4	A3070V
hCoV-19/USA/CruiseA-9/2020	EPI_ISL_413614	Nasopharyngeal swab	17 Feb	No change	NA	NA	NA
hCoV-19/USA/CruiseA-10/2020	EPI_ISL_413615	Nasopharyngeal swab	17 Feb	3259G > T	Missense	Orf1ab/NSP3	Q998H
hCoV-19/USA/CruiseA-11/2020	EPI_ISL_413616	Nasopharyngeal swab	17 Feb	10036C > T	Synonymous	Orf1ab/NSP4	NA
hCoV-19/USA/CruiseA-12/2020	EPI_ISL_413617	Oropharyngeal swab	20 Feb	29635C > T	Non-coding	3′-UTR	NA
hCoV-19/USA/CruiseA-13/2020	EPI_ISL_413618	Nasopharyngeal swab	20 Feb	No change	NA	NA	NA
hCoV-19/USA/CruiseA-14/2020	EPI_ISL_413619	Oropharyngeal swab	25 Feb	6636C > T	Missense	Orf1ab/NSP3	T2124I
				11750C > T	Missense	Orf1ab/NSP6	L3829F
				11956C > T	Synonymous	Orf1ab/NSP7	NA
hCoV-19/USA/CruiseA-15/2020	EPI_ISL_413620	Nasopharyngeal swab	18 Feb	No change	NA	NA	NA
hCoV-19/USA/CruiseA-16/2020	EPI_ISL_413621	Nasopharyngeal swab	18 Feb	No change	NA	NA	NA
hCoV-19/USA/CruiseA-17/2020	EPI_ISL_413622	Nasopharyngeal swab	24 Feb	5845A > T	Missense	Orf1ab/NSP3	K1860N
hCoV-19/USA/CruiseA-18/2020	EPI_ISL_413623	Nasopharyngeal swab	24 Feb	508_522del	Deletion	Orf1ab/NSP1	Deletion 82–86
				22033C > A	Missense	S	F157L
hCoV-19/USA/CruiseA-19/2020	EPI_ISL_414479	Nasopharyngeal swab	18 Feb	No change	NA	NA	NA
hCoV-19/USA/CruiseA-21/2020	EPI_ISL_414480	Oropharyngeal swab	21 Feb	254C > T	Non-coding	5′-UTR	NA
				29726G > T	Non-coding	3′-UTR	NA
				29751G > C	Non-coding	3′-UTR	NA
hCoV-19/USA/CruiseA-22/2020	EPI_ISL_414481	Nasopharyngeal swab	21 Feb	No change	NA	NA	NA
hCoV-19/USA/CruiseA-23/2020	EPI_ISL_414482	Nasopharyngeal swab	18 Feb	254C > T	Non-coding	5′-UTR	NA
				9157T > C	Synonymous	Orf1ab/NSP4	NA
				11083G > T	Missense	Orf1ab/NSP6	L3606F
				22104G > T	Missense	S	G181V
				29736G > T	Non-coding	3′-UTR	NA
				29751G > C	Non-coding	3′-UTR	NA
hCoV-19/USA/CruiseA-24/2020	EPI_ISL_414483	Oropharyngeal swab	17 Feb	3099C > T	Missense	Orf1ab/NSP3	T945I
				10507C > T	Synonymous	Orf1ab/NSP5	NA
				11083G > T	Missense	Orf1ab/NSP6	L3606F
				28378G > T	Synonymous	N	NA
hCoV-19/USA/CruiseA-25/2020	EPI_ISL_414484	Nasopharyngeal swab	17 Feb	11083G > T	Missense	Orf1ab/NSP6	L3606F
hCoV-19/USA/CruiseA-26/2020	EPI_ISL_414485	Oropharyngeal swab	24 Feb	11083G > T	Missense	Orf1ab/NSP6	L3606F
				28916G > A	Missense	N	G215S
hCoV-19/Japan/TK/20–31–3/2020	EPI_ISL_413459	Bronchial and lung autopsy	20 Feb	11083G > T	Missense	Orf1ab/NSP6	L3606F
				29635C > T	Non-coding	3′-UTR	NA
hCoV-19/Japan/Hu_DP_Kng_19–027/2020	EPI_ISL_412969	Throat swab	10 Feb	11083G > T	Missense	Orf1ab/NSP6	L3606F
				2963C > T	Non-coding	3′-UTR	NA
hCoV-19/Japan/Hu_DP_Kng_19–020/2020	EPI_ISL_412968	Throat swab	10 Feb	11083G > T	Missense	Orf1ab/NSP6	L3606F

del: deletion; E: envelope protein; GISAID: Global Initiative on Sharing All Influenza Data; N: nucleocapsid phosphoprotein; NA: not applicable; ORF: open reading frame; S: spike glycoprotein; SARS-CoV-2: severe acute respiratory syndrome coronavirus 2; UTR: untranslated region.

a We compared the genomes to SARS-CoV-2 sequence to WIV04 (accession no. MN996528).

Note: The genome sequences are available at Global Initiative on Sharing All Influenza Data.

Fig. 1

Phylogenetic tree of SARS-CoV-2 genomes from samples taken from people on shipboard quarantine, February 10 to February 25, 2020

SARS-CoV-2: severe acute respiratory syndrome coronavirus 2.

Notes: We generated the phylogenetic tree of the viral sequences from the cruise in MAFFT using the neighbour-joining method. Mutations appearing in more than one variant are colour coded. We have included UPHL-01 (accession no. EPI_ISL_415539U) as it is the first sample with both 11083G > T and 26326C > T.

Two mutations had previously been identified in other SARS-Cov-2 isolates. First, the deletion in hCoV-19/USA/CruiseA-18/2020 has been detected in USA-CA6 (MT044258). However, the possibility that hCoV-19/USA/CruiseA-18/2020 contributed to a new mutation in the cruise samples can be dismissed because this deletion was absent in the other samples. Second, the nonsynonymous mutation 11083G > T was found in USA-AZ1 (EPI_ISL_406223) on 22 January. Later, this mutation was also present in WA3-UW1 (EPI_ISL_413025), NY-NYUMC1 (EPI_ISL_414639) and UPHL-01 (EPI_ISL_415539U).

Viral origin and transmission

Fig. 1 shows the phylogenetic trees of the 28 samples without introducing the outgroup virus. Half (14/28) of the samples contained a virus variant that had evolved, that is, had more than two mutations, during the shipboard quarantine. We identified five subgroups after rooting the phylogenetic tree with an outgroup virus sequence of SARS-WIV16: (i) 3099C > T and 28378G > T (two samples); (ii) 11083G > T (three samples); (iii) 11410G > A and 26326C > T (two samples); (iv) 29635C > T (four samples); and (v) 29736G > T (three samples; Fig. 2). Due to 11083G > T mutation, the clustering of taxa on viral phylogenies was obvious with spatially structured host population between all these subgroups (50%; 14/28; Fig. 2).

Fig. 2

Rooted phylogenetic tree of SARS-CoV-2 genomes from samples taken from people on shipboard quarantine, February 10 to February 25, 2020

SARS-CoV-2: severe acute respiratory syndrome coronavirus 2.

Notes: Alignments of viral sequences were generated using MAFFT, and the phylogenetic tree was visualized using FigTree. Rooting was done by introducing bat SARS-like coronavirus WIV16 (accession no. KT444582) as an outgroup virus. Each colour represents a subgroup.

Whether the single 11083G > T substitution spontaneously occurred during the quarantine or the patients had been infected with a viral variant containing this mutation before boarding the ship is unclear. Nevertheless, all of the viral sequences were more similar to the WIV04 sequence than with other 143 SARS-CoV-2 isolates in the GISAID database (data repository). This result suggests that the 24 new mutations identified were generated de novo on the ship rather than deriving from multiple geographical origins.

The analysis revealed two possibilities of the viral origin: either the virus (except hCoV-19/USA/CruiseA-18/2020) originated from a single primary case with the WIV04 sequence and all substitution mutations occurred during the quarantine; or there were two simultaneously primary cases, one identical to the WIV04 sequence and one containing the 11083G > T substitution.

Natural selection of mutations

Four variants had three mutations in their genome, one variant had four mutations and two variants had six mutations (Fig. 1 and Table 1). To test the hypothesis that the virus mutation evolved under selection pressure as opposed to neutral evolution (random) onboard the cruise ship, we calculated Tajima’s D to test neutrality of DNA polymorphisms., The Tajima’s D value in the cruise was −2.03 (P < 0.01) compared with 6.75 among 143 full-length genomic sequences of SARS-CoV-2 isolates sequenced between 10 January and 13 March 2020. This result indicates that, while other SARS-CoV-2 isolates faced balancing selection (due to a positive Tajima’s D), the virus spreading on the cruise had evolved under two possible but not exclusive forces generated during the quarantine process: purifying or positive selection; and population growth of the new virus variant among infected patients after a recent bottleneck caused by the quarantine.

We further investigated these two possibilities using two neutrality tests., First, Fu and Li’s test generated a negative D value of −2.66 (P < 0.01), suggesting that the quarantine procedure provided a purifying or positive selection pressure to generate an excess of singleton sites. This conclusion was corroborated by the fact that the 39.3% (11/28) of cases contained 15 new singleton mutations (Fig. 1). Second, Zeng’s E value of  −2.37 (P < 0.01) supported the possibility of population growth of the virus after a recent bottleneck as a force. We conclude that SARS-CoV-2 viral evolution was positively correlated to the increase of the selection pressure during the shipboard quarantine.

RNA recombination

Seven samples contain the 11083G > T mutation although they belong to different subgroups (Fig. 1). By assuming the substitution rate of 0.92 × 10−3/site/year, it is unlikely that the virus variants of different subgroups all generated the same spontaneous mutation at the G11083 site into the nucleotide T within three weeks. One hypothesis is that RNA recombination occurred in these cases to gain the 11083G > T mutation.

To determine whether four variants from different subgroups (samples /USA/CruiseA-23/, /USA/CruiseA-24/, /Japan/TK/20–31–3/2020 and /Japan/Hu_DP_Kng_19–027/2020) obtained the G11083 > T mutation via RNA recombination, we analysed the patterns of linkage disequilibrium between variants with minor alleles of two SARS-CoV-2 variants (available in the data repository). Table 2 shows five pairs of mutations among all analysed samples with high log of the odds of there being a disequilibrium between two loci (> 3.0) and their r2 values. The linkages between double mutations in variants /USA/CruiseA-1/ and /USA/CruiseA-24/, /USA/CruiseA-6/ and /USA/CruiseA-7/, and /USA/CruiseA-21/ and /USA/CruiseA-23/ were statistically significant (P < 0.001). Furthermore, there was indication of positive RNA recombination. The upper 95% confidence bound of D’ for C254/G11083 was 0.86, for G29736/G11083 was 0.73, for G29751/G11083 was 0.86, for C3099/G11083 was 0.86, for G28378/G11083 was 0.86 and for C29635/G11083 was 0.80 (Fig. 3). This result supports the hypothesis that 11083G > T mutation had been gained via RNA recombination in these four variants.

Table 2

Linkage disequilibrium of SARS-CoV-2 mutations with high log of the odds and the r2 values

Pair	Log of the odds	r2
C3099/G28378	4.60	1.00
G11410/C26326	3.77	0.66
G29736/29751	3.77	0.66
C254/G29736	3.77	0.66
C254/G29751	4.60	1.00

Fig. 3

Haplotype block organization of SARS-CoV-2 mutations in samples from the cruise ship Diamond Princess, 2020

SARS-CoV-2: severe acute respiratory syndrome coronavirus 2.

Notes: Linkage disequilibrium plot of HaploView to display the confidence bounds colour scheme. Each box represents a pair of mutations. The solid spines of strong linkage disequilibrium running from one marker to another along the legs of the triangle in the linkage disequilibrium chart determine the haplotype block. We defined strong evidence of recombination if pairs for which the upper confidence bound of the coefficient of linkage disequilibrium is less than 0.9 (white squares). We did not detect any strong linkage.

We detected the mutations 11083G > T and 26326 C > T, present in the UPHL-01 sequence, in the two variants: /USA/CruiseA-6/ and /USA/CruiseA-7/. The upper 95% confidence bounds on D’ of C26326/G11083 was 0.73 (Fig. 3), suggesting that a recent RNA recombination event may have also occurred in UPHL-01.

Fig. 4 shows each haplotype in the two identified blocks and their population frequencies and connections between blocks. The value of multi-allelic D' was 0.65, which represented the level of recombination between the two blocks. These data provide evidence that RNA recombination contributed to SARS-CoV-2 mutations during shipboard quarantine.

Fig. 4

Haplotype frequencies of SARS-CoV-2 mutations in the cruise ship Diamond Princess, 2020

SARS-CoV-2: severe acute respiratory syndrome coronavirus 2.

Notes: Numbers next to each haplotype block are haplotype frequencies. The joining lines represent combined haplotypes. In the crossing areas between haplotype blocks, a value of multi-allelic D’, that is the normalized value of the coefficient of linkage disequilibrium, is shown to represent the level of recombination between blocks.

Protein mutations

The viral populations in the cruise also generated 11 new variants in the viral proteins, including proteins in the gene ORF1ab (nonstructural protein NSP3, NSP4, NSP6), and in the viral structural spike and nucleocapsid proteins (Table 1). We found five missense variants in NSP3: T945I; Q998H; P1158S in the macro domain; K1860N in papain-like protease; and T2124I in the group 2 marker domain. In the other NSPs, we detected two mutations: A3070V in NSP4 and L3829F in NSP6. For the spike protein, two mutations, F157L and G181V, were identified. For the nucleocapsid phosphoprotein, we found mutations in the RNA binding domain (P46S) and the arginine-serine domain (G215S).

Only variant /USA/CruiseA-14/ had two mutations, T2124I and L3829F in NSP3 and NSP6, respectively (Table 1).

3’-UTR mutations

In three samples, we found two mutations, 29736G > T and 29751G > T in the stem loop-II motif (Fig. 5). We used the published three-dimensional crystal structure of SARS-CoV stem loop-II motif RNA to map the nucleotides G29736 and G29751. We found that these nucleotides are equal to the nucleotides G13 and G28 in the SARS-CoV stem loop-II motif (Fig. 5). In SARS-CoV, G13 (G29736 in SARS-CoV-2) forms a base triple with A38 and C39 in a seven-nucleotide asymmetric bubble, while G28 (G29751 in SARS-CoV-2) participates in formation of an essential RNA base quartet composed of two G–C pairs (G19, C20, G28, C31; Fig. 5).

Fig. 5

Mutations in the SARS-CoV-2 stem loop-II RNA motif

SARS-CoV: severe acute respiratory syndrome coronavirus.

Notes: We compared hCoV-19/USA/CruiseA-2/2020 (accession no. EPI_ISL_413607), hCoV-19/USA/CruiseA-21/2020 (EPI_ISL_414480), hCoV-19/USA/CruiseA-23/2020 (EPI_ISL_414482) and other cruise isolates with WIV04 (MN996528) and SARS-CoV (NC004718). Conventional RNA helical base-pairing are indicated in italics. The schematic representation of the stem loop-II RNA secondary structure of SARS-CoV virus shows nucleotide interactions. The digits 13 and 28 show the location of the mutations identified in some of the cruise ship isolates and digit 19 represents the location of stem loop-II mutation found in Australia.

Discussion

Viral phylogenetics is a useful tool to study epidemiological and evolutionary processes, such as epidemic spread and spatiotemporal dynamics including metapopulation dynamics, zoonotic transmission, tissue tropism and antigenic drift. Here we report the viral phylodynamics of SARS-CoV-2 from patients in a shipboard quarantine for three weeks in February 2020. The transmission started from either one or two primary cases with WIV04 sequence and/or 11083G > T mutation, then quickly separated into at least five subgroups based on new mutations. Increased positive selection as well as RNA recombination of SARS-CoV-2 were evident during the quarantine. These results should be considered in formulation of future management protocols with respect to a SARS-CoV-2 outbreak in any relatively close quarters, such as shipboards, submarines, dormitories, prisons and hospitals.

While the quarantine averted a lot of infections on the shipboard, the phylogenetics analysis showed that viral transmission and RNA recombination occurred between the five identified subgroups. Our data fit in the coalescent model, which uses the diversity of viral genome, the viral evolutionary rate and the estimated time of infection to determine the number of viral genotypes present in the initial infected population. However, we cannot rule out that evolutionary processes, such as the transmission bottleneck that determines how much of the viral diversity generated in one host passes to another during transmission, also shaped the viral phylogenies. While spatial structure is the most general virus population structure in phylodynamic analyses, SARS-CoV-2 evolution may also have been influenced by the characteristic of the host, such as age, race and risk behaviour. Because viral transmission can preferentially occur between patients sharing any of these attributes, the real reason(s) for viral transmission between virus variants require(s) further study. Furthermore, studies on whether quarantine in close quarters also promotes virus to rapidly gain more mutations via RNA recombination are needed.

Despite the small sample size in this study, our findings from computational and statistical analyses indicated that the selection pressure was not random. We assume that SARS-CoV-2 variants were at an initial stage of evolution rather than the fixation stage, since the COVID-19 outbreak had only started 8 weeks before the Diamond Princess incident. The tools we used to measure selection pressure of adaptive evolution were Tajima’s D, Zeng’s E and Fu and Li’s D tests, because of their appropriateness in relation to the sample characteristics. For Tajima’s D, its power is not affected by RNA recombination events, number of segregation sites and various timescales nor does this test require the outgroup sequence or large sample size. Fu and Li’s D test has similar properties, but is more sensitive when singleton sites are in excess, as in our samples. As the 28 specimens were collected on the cruise ship and sequenced in Japan and the USA, the hypothesis of strict neutrality may not apply. In this case, Li et al.’s methods could improve the neutrality test for sequence polymorphisms from multiple samples. Unfortunately, we could not apply this method, since the detailed parameters needed for this method were not available in GISAID. Zeng’s E test has two advantages. First, this test is not affected by positive selection at various evolution stages but denotes powerful evidence for virus population growth. Second, the test is almost absent of any power to detect positive selection at this early stage of virus transmission, therefore the negative value is mainly contributed by the virus population growth. We excluded the method to calculate the ratio of nucleotide substitution rates, denoted dN/dS, despite the simplicity and robustness of this method. The reasons for this exclusion were: (i) the ratio is defined to describe the relative rate of selected versus neutral fixation events over long timescales, not transient polymorphisms over short timescales; and (ii) dN/dS value in a single population does not follow a monotonic function in proportion to positive selection pressure like diverged sequences.

From these analyses of selection pressure, we conclude that on the cruise, the virus evolved under strong positive selection or maybe in the process of selective sweeps, which could generate beneficial mutations for SARS-CoV-2 to quickly reach fixation.

Although RNA recombination in SARS-CoV-2 had been suggested previously, this study provides evidence that RNA recombination occurred de novo in SARS-CoV-2 genome. Within 3 weeks, the genome, sampled from four infected individuals, had gained the same 11083G > T mutation, suggesting that RNA recombination also participated in viral evolution of the virus. The RNA recombination of 11803G > T is also present in UPHL-01 variant, however, whether the carrier of the UPHL-01 acquired the variant from a cruise ship passenger or if the mutation appeared independently of cruise ship variants is unknown. Because this mutation was also later detected in other variants, in addition to UPHL-01, future studies should further investigate whether 11083G > T may increase the fitness of the carrier. Other studies have suggested the 11083G > T could be a beneficial mutation linked to asymptomatic infection.,

Of the 24 mutations we identified, 11 mutations led to amino acid substitutions and two mutations occurred in the stem loop-II motif in the 3’-UTR region. This motif is a very well conserved RNA motif in more than 30 coronaviruses., We have also reported the unique 29742G > A or 29742G > U substitutions in stem loop-II motif RNA in SARS-CoV-2 isolates in Australia (Fig. 5), reinforcing the idea that stem loop-II motif is a hotspot for mutations in SARS-CoV-2 rather than a conserved RNA domain. Mutations in this motif may disrupt RNA structure and thereby alter the viral viability or infectivity. Whether the mutations identified here, both in viral proteins and regulatory RNA regions, may enhance adaptation or attenuate virus replication or virulence requires further investigation.

We acknowledge there may be limitations in collection of patients’ samples and size in extending our findings to represent the entire SARS-CoV-2-positive population (712 people, from January 20 to March 8, 2020) in the Diamond Princess cruise. If possible, further investigation with more viral genome sequences is required to understand the detailed evolutionary lineage of transmission and whether there is difference of viral mutations between symptomatic and asymptomatic individuals.