Evolutionary history, potential intermediate animal host, and cross-species analyses of SARS-CoV-2.

To investigate the evolutionary history of the recent outbreak of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) in China, a total of 70 genomes of virus strains from China and elsewhere with sampling dates between 24 December 2019 and 3 February 2020 were analyzed. To explore the potential intermediate animal host of the SARS-CoV-2 virus, we reanalyzed virome data sets from pangolins and representative SARS-related coronaviruses isolates from bats, with particular attention paid to the spike glycoprotein gene. We performed phylogenetic, split network, transmission network, likelihood-mapping, and comparative analyses of the genomes. Based on Bayesian time-scaled phylogenetic analysis using the tip-dating method, we estimated the time to the most recent common ancestor and evolutionary rate of SARS-CoV-2, which ranged from 22 to 24 November 2019 and 1.19 to 1.31 × 10-3 substitutions per site per year, respectively. Our results also revealed that the BetaCoV/bat/Yunnan/RaTG13/2013 virus was more similar to the SARS-CoV-2 virus than the coronavirus obtained from the two pangolin samples (SRR10168377 and SRR10168378). We also identified a unique peptide (PRRA) insertion in the human SARS-CoV-2 virus, which may be involved in the proteolytic cleavage of the spike protein by cellular proteases, and thus could impact host range and transmissibility. Interestingly, the coronavirus carried by pangolins did not have the RRAR motif. Therefore, we concluded that the human SARS-CoV-2 virus, which is responsible for the recent outbreak of COVID-19, did not come directly from pangolins.

INTRODUCTION

On 11 February 2020, the International Committee on Taxonomy of Viruses officially renamed the novel coronavirus (ie, previously 2019‐nCoV) responsible for the current outbreak of COVID‐19, severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2). This was chosen based on analysis of the new coronavirus's evolutionary history and associated pathogen (ie, SARS‐CoV). The virus, which emerged in December 2019 in the Chinese city of Wuhan, causes a respiratory illness called COVID‐19, which can spread from person to person. , As of 21 February 2020, there have been 76 288 cases of SARS‐CoV‐2 confirmed in mainland China, including 11 477 serious, 2 345 deaths, and 20 659 discharged, as well as 68 cases in Hong Kong, 10 in Macao, and 26 in Taiwan. More than 1300 cases have also been confirmed in at least 27 other countries on four continents. World Health Organization (WHO) officials outlined their top research priorities for controlling the outbreak of the coronavirus‐associated disease known as COVID‐19 and highlighted the importance of developing candidate therapeutics and easy‐to‐apply diagnostics for identifying active, asymptomatic, and resolved infections. Of note, the Coronaviridae family not only includes SARS‐CoV‐2, but also SARS‐CoV, Middle East respiratory syndrome coronavirus (MERS‐CoV), and common cold viruses (eg, 229E, OC43, NL63, and HKU1). The SARS‐CoV pathogen was responsible for >8 000 cases and 774 deaths in 37 countries during the 2002 to 2003 SARS outbreak, , , and the MERS‐CoV pathogen was responsible for 2 494 cases and 858 deaths in 27 countries during the 2012 MERS outbreak. ,

Coronaviruses are known to circulate in mammals and birds. Previous studies revealed that both SARS‐CoV and MERS‐CoV are zoonotic in origin, originally coming from bats, , , , with SARS‐CoV spreading from bats to palm civets to humans, , , and MERS‐CoV spreading from bats to camels to humans. , Recent research has also reported that the SARS‐CoV‐2 virus likely originated in bats, a proposal based on the similarity of its genetic sequence to that of other known coronaviruses. However, like SARS‐CoV, MERS‐CoV, and many other coronaviruses, the SARS‐CoV‐2 virus may have been transmitted to humans by an intermediate animal host. Therefore, the identity of the animal source of SARS‐CoV‐2 remains a key and urgent question. Furthermore, to stem future outbreaks of this type and preventing the transmission of zoonotic diseases to humans should be a top research priority.

The existence of an intermediate animal host of SARS‐CoV‐2 between a probable bat reservoir and humans is still under investigation. The discovery of a virus closely related to the newly emerged SARS‐CoV‐2 in a data set from pangolins sampled more than a year ago illustrates that the sampling of other mammals handled or consumed by humans could uncover even more closely related viruses.

During a press conference on 7 February 2020, two researchers (Shen Yongyi and Xiao Lihua) from the South China Agricultural University in Guangzhou identified the pangolin as a potential source of the SARS‐CoV‐2 virus based on genetic comparison of coronaviruses taken from pangolins and from humans infected during the recent outbreak. By analyzing more than 1 000 metagenomic samples and using molecular biology testing, they found that the positive rate of β coronavirus in pangolins was 70% and that the genome sequence of an isolated virus strain was 99% similar to that of the SARS‐CoV‐2 virus. Thus, whether pangolins acted as a direct intermediate animal host of the SARS‐CoV‐2 virus is worth further investigation.

In the present study, we performed analyses of the transmission dynamics and evolutionary history of the virus based on 70 genomes of SARS‐CoV‐2 strains sampled from Australia (n = 4), Belgium (n = 1), China (Hubei Province, n = 19; Guangdong Province, n = 16; Zhejiang Province, n = 4; Taiwan, n = 1), Finland (n = 1), France (n = 4), Germany (n = 1), Japan (n = 1), Korea (n = 1), Singapore (n = 3), Thailand (n = 2), United Kingdom (n = 2), and United States (n = 10) with sampling dates between 24 December 2019 and 3 February 2020. We reanalyzed two of the 21 pangolin metagenome samples from previously published data and compared the amino acid sequences of the S protein of SARS‐CoV‐2 and SARS‐related coronaviruses. These analyses should extend our understanding of the origins and dynamics, cross‐species transmission, and subsequent host adaptation of the SARS‐CoV‐2 outbreak in China and elsewhere.

MATERIALS AND METHODS

Collation of SARS‐CoV‐2 genome data sets

As of 9 February 2020, 73 genomes of SARS‐CoV‐2 strains obtained from humans have been released on GISAID (http://gisaid.org/). The BetaCoV/Wuhan/IPBCAMS‐WH‐02/2019 (EPI_ISL_403931), BetaCoV/Shenzhen/SZTH‐001/2020 (EPI_ISL_406592), and BetaCoV/Shenzhen/SZTH‐004/2020 (EPI_ISL_406595) samples show evidence of sequencing artefacts due to the appearance of clustered spurious single nucleotide polymorphisms (SNPs) and were thus excluded from this study. The final data set (“dataset_70”) included 70 genomes of SARS‐CoV‐2 strains from Australia (n = 4), Belgium (n = 1), China (n = 40), Finland (n = 1), France (n = 4), Germany (n = 1), Japan (n = 1), Korea (n = 1), Singapore (n = 3), Thailand (n = 2), UK (n = 2), and USA (n = 10) with sampling dates between 24 December 2019 and 3 February 2020. Of the 40 samples collected from China, 19 were from Hubei Province, 16 were from Guangdong Province, 4 were from Zhejiang Province, and 1 was from Taiwan (Table S1).

To investigate the potential intermediate hosts of SARS‐CoV‐2 (between originating animal and human hosts), two samples (SRR10168377 and SRR10168378) obtained from previously reported Malayan pangolin (Manis javanica) viral metagenomic sequencing data (Bio Project PRJNA573298) were downloaded from the NCBI SRA public database. After assembly, the SRR10168377 and SRR10168378 genomes were 16 999 and 6 392 bp in length, respectively. We defined another data set (“dataset_6”) composed of six genome sequences of coronavirus strains. BetaCoV/Wuhan‐Hu‐1/2019 (EPI_ISL_402125) was grouped as “Clade A,” one (BetaCoV/bat/Yunnan/RaTG13/2013; EPI_ISL_402131) and two (bat‐SL‐CoVZC45; MG772933 and bat‐SL‐CoVZXC21; MG772934) SARS‐related coronaviruses were grouped as “Clade B” and “Clade D,” respectively. The two assembled genomes from SRR10168377 and SRR10168377 were grouped into “Clade C.” The two data sets (“dataset_70” and “dataset_6”) were aligned using MAFFT v7.222 and then manually curated using BioEdit v7.2.5.

Recombination and phylogenetic analyses

To assess the recombination of “dataset_70,” we employed the pairwise homoplasy index (PHI) to measure the similarity between closely linked sites using SplitsTree v4.15.1. The best‐fit nucleotide substitution models for the two data sets were identified according to the Bayesian information criterion (BIC) method with 3 (24 candidate models) or 11 (88 candidate models) substitution schemes in jModelTest v2.1.10.  To evaluate the phylogenetic signals of “dataset_70” and “dataset_6,” we performed likelihood‐mapping analysis using TREE‐PUZZLE v5.3, with 25 000 to 175 000 randomly chosen quartets for the two data sets. For “dataset_70,” split network analysis was performed using Kishino‐Yano‐85 distance transformation with the NeighborNet method, which can be loosely thought of as a “hybrid” between the neighbor‐joining (NJ) and split decomposition methods, implemented in TREE‐PUZZLE v5.3. For “dataset_70,” NJ phylogenetic trees were constructed using the Kimura 2‐parameter method implemented in MEGA v7.0.26. For “dataset_6,” NJ phylogenetic trees were constructed using the Maximum Composite likelihood (MCL) method, and rate variation among sites was modeled with a gamma distribution (shape parameter = 4) in MEGA v7.0.26. For “dataset_70,” maximum‐likelihood (ML) phylogenies were reconstructed using the Hasegawa‐Kishino‐Yano (HKY) nucleotide substitution model in PhyML v3.1. For “dataset_6,” ML phylogenies were reconstructed using the general time reversible nucleotide substitution model with gamma‐distributed rate variation among sites (GTR + G) model in PhyML v3.1. For all NJ and ML phylogenies of the two data sets, bootstrap support values were calculated with 1 000 replicates and trees were midpoint rooted. For “dataset_70,” regression analyses were used to determine the correlations among sampling dates and root‐to‐tip genetic divergences of the respective ML phylogenies with TempEst v1.5. We also estimated the evolutionary rate and time to the most recent common ancestor (TMRCA) for “dataset_70” using ML dating in the TreeTime package.

Reconstruction of time‐scaled phylogenies

To reconstruct the evolutionary history of SARS‐CoV‐2, Bayesian inference through a Markov chain Monte Carlo (MCMC) framework was implemented in BEAST v1.8.4, with the BEAGLE v2.1.2 library program used for computational enhancement. We used two schemes to set the time‐scale prior for each data set: that is, constrained evolutionary rate method with a log‐normal prior (mean = 1.0 × 10−3 substitutions per site per year; 95% Bayesian credible interval [BCI]: 1.854 × 10−4 to 4 × 10−3 substitutions per site per year) placed on the evolutionary rate parameter, as per previous studies, , , and the tip‐dating method, for which the overall estimated evolutionary rate was given an uninformative continuous‐time Markov chain (CTMC) reference prior. We ran Bayesian phylogenetic analyses using various clock model combinations (ie, strict clock and uncorrelated log‐normal relaxed clock ) and coalescent tree priors (ie, constant size and exponential growth). To ensure adequate mixing of model parameters, MCMC chains were run for 100 million steps with sampling every 10 000 steps from the posterior distribution. Convergence was evaluated by calculating the effective sample sizes of the parameters using Tracer v1.7.1. All parameters had an effective sample size >200, indicative of sufficient sampling. Trees were summarized as maximum‐clade credibility (MCC) trees using TreeAnnotator v1.8.4 after discarding the first 10% as burn‐in, and then visualized in FigTree v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree).

Transmission network reconstruction

The HIV TRAnsmission Cluster Engine (HIV‐TRACE; www.hivtrace.org) was employed to infer transmission network clusters for SARS‐CoV‐2 “dataset_70.” All pairwise distances were calculated and the putative linkages between each pair of genomes were considered whenever their divergence was ≤0.0001 (0.01%) or ≤0.00001 (0.001%) substitutions/site (TN93 substitution model). Multiple linkages were then combined into putative transmission clusters. Clusters that comprised of only two linked nodes were identified as dyads. This approach detects transmission clusters in which the clustering strains are genetically similar, implying a direct or indirect epidemiological connection.

Similarity plot analysis

To investigate the putative parents of SARS‐CoV‐2, we performed similarity plot analysis based on the Kimura two‐parameter method with a window size of 200 bp and step size of 20 bp using SimPlot v.3.5.14. We divided “dataset_6” into four clades (ie, Clade A, Clade B, Clade C, and Clade D), with Clade A designated as the query group.

RESULTS

Demographic characteristics of SARS‐CoV‐2

“Dataset_70” included 70 genomes of SARS‐CoV‐2 strains sampled from Australia (n = 4), Belgium (n = 1), China (Hubei Province, n = 19; Guangdong Province, n = 16; Zhejiang Province, n = 4; Taiwan, n = 1), Finland (n = 1), France (n = 4), Germany (n = 1), Japan (n = 1), Korea (n = 1), Singapore (n = 3), Thailand (n = 2), UK (n = 2), and USA (n = 10) with sampling dates between 24 December 2019 and 3 February 2020 (Table S1). The samples were primarily from China (57.14%) and Hubei Province (27.14%), the Chinese Province was acknowledged as the original epicenter of the SARS‐CoV‐2 outbreak.

Tree‐like signals and phylogenetic analyses

For “dataset_70” and “dataset_6,” HKY and GTR + G were the models of best‐fit, respectively, across the two different substitution schemes (ie, 24 and 88 candidate models) according to the BIC method, and were thus used in subsequent likelihood‐mapping and phylogenetic analyses for the two data sets. The PHI tests of “dataset_70” did not find statistically significant evidence of recombination (P = 1.0). Likelihood‐mapping analysis of “dataset_70” revealed that 69.7% of the quartets were distributed in the center of the triangle, indicating a strong star‐like topology signal reflecting a novel virus, which may be due to exponential epidemic spread (Figure 1A). Likewise, 25.9% of the quartets from “dataset_6” were distributed in the center of the triangle, indicating a strong phylogenetic signal (Figure 1B). The split network generated for “dataset_70” using the NeighborNet method was highly unresolved, and the phylogenetic relationship of “dataset_70” was probably best represented by a network rather than a tree (Figure 1C). The existence of polytomies indicated—in contrast to that expected in a strictly bifurcating tree—an explosive, star‐like evolution of SARS‐CoV‐2. Both the NJ and ML phylogenetic analyses of SARS‐CoV‐2 “dataset_70” also showed star‐like topologies, in accordance with the likelihood‐mapping results (Figures 2 and S1). The ML phylogenetic tree showed greater star‐like topology than the NJ phylogenetic tree, indicating that the ML method was more reasonable for “dataset_70.” Root‐to‐tip regression analyses between genetic divergence and sampling date using the best‐fitting root showed that “dataset_70” had a minor strong positive temporal signal (R  = 0.0808; correlation coefficient = 0.2843) (Figure 3). This result suggests a minor clock‐like pattern of molecular evolution, with an estimated substitution rate of 3.3452 × 10−4 substitutions per site per year and TMRCA occurring on 19 October 2019. ML dating analyses between genetic divergence and sampling date also showed that “dataset_70” had a minor strong positive temporal signal (R = 0.08) (Figure S2). The estimated evolutionary rate and TMRCA were 3.34 × 10−4 substitutions per site per year and 19 October 2019, respectively, in accordance with the root‐to‐tip regression results. Based on Bayesian time‐scaled phylogenetic analysis using the constrained evolutionary rate method with a log‐normal prior (mean = 1.0 × 10−3 substitutions per site per year; 95% BCI: 1.854 × 10−4 to 4 × 10−3 substitutions per site per year) placed on the evolutionary rate parameter, the estimated TMRCA dates and evolutionary rates for SARS‐CoV‐2 from “dataset_70” ranged from 21 May 2019 to 13 October 2019 (95% BCI: 27 and 30 January 2020) and from 1.57 × 10−4 to 1.06 × 10−3 substitutions per site per year (95% BCI: 1.08 × 10−4 to 3.10 × 10−3), respectively (Table 1). Furthermore, based on Bayesian time‐scaled phylogenetic analysis using the tip‐dating method, the estimated TMRCA dates and evolutionary rates from “dataset_70” ranged from 22 to 24 November 2019 (95% BCI: 23 October 2019 and 16 December 2019) and from 1.19 × 10−3 to 1.31 × 10−3 substitutions per site per year (95% BCI: 6.22 × 10−4 to 1.96 × 10−3), respectively (Table 1). Thus, the estimated TMRCA dates and evolutionary rates for SARS‐CoV‐2 from “dataset_70” were consistent among the different clock models (strict and relaxed) but were distinct among the different dating methods (constrained‐dating and tip‐dating). The estimated TMRCA dates and evolutionary rates for SARS‐CoV‐2 from “dataset_70” using the tip‐dating method exhibited much narrower 95% BCIs than the constrained‐dating method. In addition, the estimated TMRCA dates and evolutionary rates for SARS‐CoV‐2 from “dataset_70” were consistent between the different coalescent tree models (ie, constant and exponential) when using the tip‐dating method but were distinct when using the constrained‐dating method. For each data set, we employed the HKY nucleotide substitution model, as well as a constant size coalescent tree prior and strict molecular clock model to estimate the TMRCA. The estimates of the MCC phylogenetic relationships among the SARS‐CoV‐2 genomes from the Bayesian coalescent framework using the tip‐dating method, as well as the constant size coalescent tree prior and strict molecular clock, are displayed in Figure 4. As shown, eight phylogenetic clusters (number of sequences 2 to 7; posterior probability 0.99 to 1.0) were identified.

Figure 1

Likelihood‐mapping and split network analyses of severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2). Likelihoods of three tree topologies for each possible quartet (or for a random sample of quartets) are denoted by data points in an equilateral triangle. Distribution of points in seven areas of the triangle reflects tree‐likeness of data. Specifically, three corners represent fully resolved tree topologies; center represents an unresolved (star) phylogeny; and sides represent support for conflicting tree topologies. Results of likelihood‐mapping analyses of two data sets (“dataset_70,” A; and “dataset_6,” B) and split network analyses of “dataset_70” (C) are shown

Figure 2

Estimated maximum‐likelihood phylogenies of severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2). Colors indicate different sampling locations. Tree is midpoint rooted. Results of maximum‐likelihood phylogenetic analyses of “dataset_70” are shown

Figure 3

Regression of root‐to‐tip genetic distance against year of sampling for severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2). Colors indicate different sampling locations. Gray indicates linear regression line. Results of linear regression analyses of “dataset_70” are shown

Table 1

Bayesian phylogenetic estimates of evolutionary parameters for genome sequences of 2019‐novel coronavirus (2019‐nCoV) under different clock models and coalescent tree priors

Clock model	Clock prior	Coalescent tree prior	Substitution rate (substitutions per site per year)			Clade A MRCA
			Mean	Lower 95% HPD	Upper 95% HPD	Mean	Lower 95% HPD	Upper 95% HPD
Strict	Constraint‐dating	Constant	1.06E‐03	1.24E‐04	3.10E‐03	2019‐10‐13	2019‐04‐05	2020‐01‐30
		Exponential	1.58E‐04	1.10E‐04	2.27E‐04	2019‐05‐23	2019‐01‐31	2019‐09‐06
	Tip‐dating	Constant	1.24E‐03	6.74E‐04	1.82E‐03	2019‐11‐24	2019‐10‐29	2019‐12‐15
		Exponential	1.19E‐03	6.22E‐g04	1.81E‐03	2019‐11‐22	2019‐10‐23	2019‐12‐15
Relaxed	Constraint‐dating	Constant	1.01E‐03	1.11E‐04	3.01E‐03	2019‐09‐29	2019‐03‐05	2020‐01‐29
		Exponential	1.57E‐04	1.08E‐04	2.26E‐04	2019‐05‐21	2019‐01‐27	2019‐09‐05
	Tip‐dating	Constant	1.31E‐03	7.40E‐04	1.96E‐03	2019‐11‐24	2019‐10‐25	2019‐12‐16
		Exponential	1.28E‐03	6.46E‐04	1.92E‐03	2019‐11‐23	2019‐10‐24	2019‐12‐16

Abbreviations: 2019‐nCoV, 2019‐novel coronavirus; HPD, highest posterior density; MRCA, most recent common ancestor.

Figure 4

Estimated maximum‐clade‐credibility tree of severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2) using tip‐dating method. Colors indicate different sampling locations. Nodes are labeled with posterior probability values. Estimated maximum‐clade credibility (MCC) tree of “dataset_70” are shown

Transmission network analysis

We considered individuals as genetically linked when the genetic distance between SARS‐CoV‐2 strains was <0.01% substitutions/site. Based on this, we identified one large transmission cluster that included 66 of 70 (94.29%) genomes, thus suggesting low genetic divergence for “dataset_70” (Figure S3). We also considered individuals as genetically linked when the genetic distance between SARS‐CoV‐2 strains was <0.001% substitutions/site. Based on this, we identified 6 transmission clusters that included 37 of 70 (52.86%) genomes for “dataset_70” (Figure 5). Clusters ranged in size from 2 to 23 genomes.

Figure 5

Transmission clusters of severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2). Structure of inferred SARS‐CoV‐2 transmission clusters from “dataset_70” using genetic distances of <0.001% substitutions/site is shown. Nodes (circles) represent connected individuals in overall network, and putative transmission linkages are represented by edges (lines). Nodes are color‐coded by sampling locations

Potential intermediate host analyses for SARS‐CoV‐2

The NJ and ML phylogenetic topologies of “dataset_6” were consistent with each other (Figure S4), indicating that the use of a small number of sequences could show similar topological results. Homology plot analysis of “dataset_6” also revealed that BetaCoV/bat/Yunnan/RaTG13/2013 was more similar to the SARS‐CoV‐2 virus than the coronavirus obtained from the two pangolin samples (SRR10168377 and SRR10168378), consistent with phylogenetic analysis (Figure S5). Of note, “Clade D” (bat‐SL‐CoVZC45 and bat‐SL‐CoVZXC21) had higher similarity to the SARS‐CoV‐2 virus in the first 12 000 bp region of the full alignment than to the pangolin coronavirus (Figure S5). We also found that a unique peptide (PRRA) insertion region in the spike protein at the junction of S1 and S2 junction in the human SARS‐CoV‐2 virus (“Clade A”) induced a furin cleavage motif (RRAR), which could be a predicted polybasic cleavage site, and thus a unique feature of SARS‐CoV‐2, in comparison with the other three clades (“Clade B,” “Clade C,” and “Clade D”) (Figure S6).

DISCUSSION

On the basis of “dataset_70,” our likelihood‐mapping analysis confirmed additional tree‐like signals over time compared with our previous results. , This result implies increasing genetic divergence of SARS‐CoV‐2 in human hosts (Figure 1A), consistent with the findings of our earlier studies. , Split network analysis for SARS‐CoV‐2 “dataset_70” using the NeighborNet method was highly unresolved, indicating an explosive, star‐like evolution of SARS‐CoV‐2, and recent and rapid human‐to‐human transmission (Figure 1C). These results are consistent with the ML phylogenetic analyses, which showed polytomy topology from “dataset_70” (Figure 2). However, NJ phylogenetic analyses showed a more bifurcating tree topology compared with the ML phylogenetic analyses (Figure S1). This is a good example showing the differences between NJ and ML phylogenetic construction methods. “Dataset_70” had a minor strong positive temporal signal based on root‐to‐tip regression and ML dating analyses (Figures 3 and S2), with the estimated TMRCA dates and evolutionary rates for SARS‐CoV‐2 found to be nearly identical using both analyses (Table 1). The estimated TMRCA dates and evolutionary rates for SARS‐CoV‐2 were very similar across different clock models and coalescent tree priors using the tip‐dating method. The estimated TMRCA dates and evolutionary rates for SARS‐CoV‐2 were also very similar across different clock models using the constrained‐dating method, but highly distinct across the different coalescent tree priors (Table 1). The TMRCA estimated by the tip‐dating method was relatively narrower than that determined by the constrained‐dating method, consistent with our previous studies. , Bayesian analyses with the tip‐dating method using a strict clock as well as constant size coalescent tree prior indicated that SARS‐CoV‐2 is evolving at a rate of 1.24 × 10−3 substitutions per site per year (Table 1), in accordance with our prior research , and similar to that found for other human coronaviruses. Our results also suggest that the virus originated on 24 November 2019, which is in further agreement with our earlier studies. , We identified eight phylogenetic clusters (number of sequences 2 to 7) with posterior probabilities between 0.99 and 1.0 using Bayesian inference (Figure 4). We also identified six transmission clusters (number of sequences 2 to 23) when the genetic distance between the SARS‐CoV‐2 strains was <0.001% substitutions/site (Figure 5). However, our conclusions should be considered preliminary and explained with caution due to the limited number of SARS‐CoV‐2 genome sequences presented in this study. As more genome sequences become available, there may be stronger among‐lineage rate variation over time as to warrant using a relaxed clock model, but we anticipate that the evolutionary rates and TMRCA dates will be broadly similar to those estimated here. As the number of substitutions is still small, it is tempting to speculate that sequencing errors could have a considerable impact on the evolutionary rate and TMRCA date estimates. We removed three SARS‐CoV‐2 genome sequences (ie, BetaCoV/Wuhan/IPBCAMS‐WH‐02/2019, EPI_ISL_403931; BetaCoV/Shenzhen/SZTH‐001/2020, EPI_ISL_406592; BetaCoV/Shenzhen/SZTH‐004/2020, EPI_ISL_406595) with potential sequencing errors, but these may have less impact on the above estimates when more substitutions of SARS‐CoV‐2 are accumulated over time. We also expect that as more SARS‐CoV‐2 genome sequences become available, the estimated 95% BCIs of the evolutionary rates and TMRCA dates will be narrower.

We found that the Pangolin‐CoV virus from the two pangolin samples was clustered with the SARS‐CoV‐2 virus with 100% bootstrap support; however, BetaCoV/bat/Yunnan/RaTG13/2013 was more similar to the SARS‐CoV‐2 virus than to the pangolin coronavirus and the human SARS‐CoV‐2 virus (“Clade A”) showed a unique peptide (PRRA) insertion not found in the other three clades (“Clade B,” “Clade C,” “Clade D”). This insertion constitutes an RRAR motif in the spike protein at the junction of S1 and S2 junction in the human SARS‐CoV‐2 virus, after considering the next amino acid (R) of the unique peptide (PRRA) (Figure S6). Of note, the highly favored motifs for furin cleavage are Arg‐X‐(Arg/Lys)‐Arg (RXRR or RXKR), and the minimal motifs for furin cleavage can be RXXR. We also note that some of the other coronaviruses have a furin motif in almost the same location in their spike proteins. , Lentiviruses have an RKXR (R, arginine; K, lysine; X, any amino acid) site between gp120 and gp41, cleaved by furin to convert gp160 into subunits. , , , Therefore, it is tempting to speculate that cleavage or lack of cleavage of the spike protein at this site could significantly impact host range and transmissibility. Taken together, the pangolin coronavirus samples (SRR10168377 and SRR10168378) were less similar to the SARS‐CoV‐2 virus than to the BetaCoV/bat/Yunnan/RaTG13/2013 virus and did not have the RRAR motif. Therefore, we concluded that the human SARS‐CoV‐2 virus, which is responsible for the current outbreak of COVID‐19, did not come directly from pangolins. However, due to the limited viral metagenomic data obtained from pangolins, we cannot exclude that other pangolins from China may contain coronaviruses that exhibit greater similarity to the SARS‐CoV‐2 virus.

In conclusion, our results emphasize the importance of further epidemiological investigations, genomic data surveillance, and experimental studies of the role of the unique furin cleavage motif (RRAR) of SARS‐CoV‐2 in the spike protein at the junctions of S1 and S2. Such work could positively impact public health in terms of guiding prevention efforts to reduce SARS‐CoV‐2 transmission in real time, and to stem future outbreaks of zoonotic diseases.

CONFLICT OF INTERESTS

The authors declare that there are no conflict of interests.

AUTHOR CONTRIBUTIONS

XL conceived and designed the study and drafted the manuscript. XL, BF, and AC analyzed the data. XL, QN, JZ, QZ, YL, BF, and AC interpreted the data and provided critical comments. All authors reviewed and approved the final manuscript.

Supporting information

Click here for additional data file.

Supporting information

Click here for additional data file.

Supporting information

Click here for additional data file.

Supporting information

Click here for additional data file.

Supporting information

Click here for additional data file.

Supporting information

Click here for additional data file.

Supporting information

Click here for additional data file.

Supporting information

Click here for additional data file.