Molecular evolutionary characteristics of SARS-CoV-2 emerging in the United States.

SARS-CoV-2 is a newly discovered beta coronavirus at the end of 2019, which is highly pathogenic and poses a serious threat to human health. In this paper, 1875 SARS-CoV-2 whole genome sequences and the sequence coding spike protein (S gene) sampled from the United States were used for bioinformatics analysis to study the molecular evolutionary characteristics of its genome and spike protein. The MCMC method was used to calculate the evolution rate of the whole genome sequence and the nucleotide mutation rate of the S gene. The results showed that the nucleotide mutation rate of the whole genome was 6.677 × 10-4 substitution per site per year, and the nucleotide mutation rate of the S gene was 8.066 × 10-4 substitution per site per year, which was at a medium level compared with other RNA viruses. Our findings confirmed the scientific hypothesis that the rate of evolution of the virus gradually decreases over time. We also found 13 statistically significant positive selection sites in the SARS-CoV-2 genome. In addition, the results showed that there were 101 nonsynonymous mutation sites in the amino acid sequence of S protein, including seven putative harmful mutation sites. This paper has preliminarily clarified the evolutionary characteristics of SARS-CoV-2 in the United States, providing a scientific basis for future surveillance and prevention of virus variants.

INTRODUCTION

Coronaviruses (CoV) is an enveloped RNA virus, widely distributed in humans and other mammals, and can cause a variety of diseases, e.g., respiratory intestinal, liver, and nervous system diseases. , , , , , There have been seven types of coronaviruses that infect humans, including SARS‐CoV, MERS‐CoV, HCoV‐229E, HCoV‐HKU1, HCoV‐NL63, HCoV‐OC43, and SARS‐CoV‐2, and the genome sequence and spike protein structure of SARS‐CoV‐2 are similar to that of SARS‐CoV. , , The receptor binding and membrane fusion are the initial and critical steps in the SARS‐CoV‐2 infection cycle, during which the SARS‐CoV‐2 spike protein (S protein) plays a key role. Therefore, it is very important to study the SARS‐CoV‐2 S protein. ,

As of July 2021, the number of SARS‐CoV‐2 infections has reached 180 million people worldwide. As for the United States, more than 33 million COVID‐19 cases have been diagnosed, which is close to a quarter of the total number of confirmed cases worldwide. The death toll in the United States exceeds 60 million.

In previous studies, Li et al. analyzed the global evolution rate of SARS‐CoV‐2 in the first month of the outbreak, with an estimated mean nucleotide mutation rate ranging from 1.7926 × 10−3 to 1.8266 × 10−3 substitution per site per year. Four months after the outbreak, the mutation rate became 3.95 × 10−4 per nucleotide per year, which was almost seven times lower than the mutation rate of SARS‐CoV and two times lower than that of MERS‐CoV. , , Motayo et al. reported that the evolution rate of SARS‐CoV‐2 in Africa from February 24 to April 24 was 4.133 × 10−4 substitution per nucleotide per year. The Nextstrain website estimates that the annual nucleotide evolution rate is 8 × 10−4 substitution per nucleotide per year based on current statistics. As for the evolution rate of the gene encoding the S protein, Pereson et al. reported that the global evolution rate was 2.19 × 10−3 nucleotide substitutions per site per year as of April 2020, while the evolution rate reduced to 1.08 × 10−3 as of September 2020. However, so far, no relevant research has been reported on the US cases, so we initiate this study.

MATERIALS AND METHODS

Sequence data collection

Various coronavirus genome sequence data for the phylogenetic analysis

A total of 15 whole genome sequences and the corresponding S gene sequences were employed in this study from the NCBI database (https://www.ncbi.nlm.nih.gov/) (Table S1), including SARS‐CoV‐2, SARS‐CoV, MERS‐CoV, HCoV‐229E, HCoV‐OC43, HCoV‐NL63, HCoV‐HKU1, Bat SARS‐like, Civet SARS‐CoV, Pangolin CoV, Murine hepatitis virus, Rat CoV, Erinaceus CoV, Camel MERS‐CoV, and Bovine CoV.

SARS‐CoV‐2 genome sequence and the corresponding S gene sequence data collected from US cases to study molecular evolution

A total of 2241 whole‐genome sequences of SARS‐CoV‐2 and the corresponding S gene sequences (from 2020/2/27 to 2021/4/8) were obtained from the severe acute respiratory syndrome coronavirus 2 data hub of NCBI Virus (https://www.ncbi.nlm.nih.gov/labs/virus/vssi/#/) (Table S2). After removing the low‐quality sequences (with more than 10 ambiguous nucleotides), 1875 SARS‐CoV‐2 genome sequences were retained for further analysis, using MN908947 (isolate from Wuhan‐Hu‐1) as a reference sequence.

In addition, another 69 098 S protein amino acid sequences obtained from the NCBI database were used for the mutation analysis of amino acids, which were also retrieved from the SARS‐CoV‐2 cases in the United States.

Sequence alignment and phylodynamics analysis

A multiple sequence alignment of the 1875 SARS‐CoV‐2 genome sequences was performed using MAFFT v7.464. ModelFinder was used to analyze the optimal substitution model for the genome sequences based on the results of Bayesian Information Criterion, with a result of the best nucleotide substitution model to be GTR+F+G4.

Mutation rate can be defined as the number of mutations per cell division, per generation or per unit time. Markov Chain Monte Carlo (MCMC) method was used to reconstruct the maximum clade credibility (MCC) tree and calculate the mutation rate by BEAST V2.6.2. To set the time scale prior for the dataset, we used a constrained evolution rate with a Log‐normal prior averaged at 10−3 by substitution per site per year. We performed a phylogenetic Bayesian analysis using the Relaxed Clock Log‐Normal molecular clock model and selected the Coalescent Bayesian Skyline as the model of population size and growth according to the relevant studies. , , The whole‐genome sequence and the sequence coding spike protein sequence (S gene) were separately analyzed by MCMC to calculate the mutation rate, with a length of 4 × 108 steps, sampling every 4 × 104 steps. The convergence of all the parameters (ESS >200, burn‐in 10%) was verified with Tracer v1.7.1. The final MCC tree was generated by TreeAnnotator (a software package in BEAST2) and displayed in Figtree v1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/).

Adaptive evolution analysis and recombination analysis

The SARS‐CoV‐2 MN908947 was used as the reference sequence, which was aligned against the genome sequences of the other types of coronavirus by MAFFT. IQ‐Tree 2.0.3 was used to establish the ML tree. The CODEML program in PAML was used for the selection pressure analysis. In the branch‐site model, the SARS‐CoV‐2 was set as the foreground and the other viruses as the background. Full genomes and S gene sequences were analyzed separately to detect recombination events using RDP4 and SimPlot. ,

The S protein nonsynonymous or indel variant biological function analysis

PROVEAN was used to predict whether nonsynonymous mutations of S protein would affect its function, which can detect harmful substitution of amino acid and predict whether the substitution will affect its phenotype. A PROVEAN score of −2.5 or lower indicates that amino acid substitution is harmful, and a score higher than −2.5 is considered neutral.

RESULTS AND DISCUSSION

Mutation rate in SARS‐CoV‐2

According to the results of BEAST2, we obtained the mutation rate of 6.677 × 10−4 per site per year (95% highest posterior density [HPD]: 6.117 × 10−4 to 7.270 × 10−4) for the whole genome and of 8.066 × 10−4 per site per year (95% HPD: 5.969 × 10−4 to 1.038 × 10−3) for the sequence coding S protein. The first more infectious mutation discovered by the scientists is D614G, which is caused by the change of base A to G at position 23403. We identified 12 high‐frequency mutations in the S gene among the 69 098 amino acid sequences, with the highest mutation frequency to be 98.47% at the position of amino acid (AA) 614 (Figure 1).

Figure 1

Distribution of mutations in the S protein. All mutations in the S protein are nonsynonymous, with position AA614 having the highest frequency

Figure 2 shows the MCC tree with Bayesian phylogeographic reconstruction of SARS‐CoV‐2 isolates. To make the image visual effect clearer without losing the representativeness of the genomes, we only show 380 genome sequences in Figure 2 to display the MCC tree, and the MCC tree including the complete list of 1875 sequences is shown in Figure S1. The detailed information obtained from BEAST2 results is shown in Table S3. We estimate that November 5, 2019, is the time of the most recent common ancestor of SARS‐CoV‐2 emerging in the United States, with the 95% HPD ranging from September 21, 2019, to December 16, 2019. This conclusion is also consistent with the first case reported in the literature on December 1, 2019.

Figure 2

Bayesian maximum clade credibility tree

Phylogenetic analysis of SARS‐CoV‐2

To explore the evolutionary selection pressure of SARS‐CoV‐2, a phylogenetic analysis was carried out using MN908947 as the reference sequence, and the whole genome of MN908947 was aligned with the other 14 coronaviruses to establish the phylogenetic tree. MN908947 was used as the foreground branch with the other 14 coronaviruses as the background branch for the branch‐site model. The analysis results are shown in Table 1 and Table S4. Thirteen and two statistically significant positive sites were detected in the whole genome and S gene respectively (p < 0.005).

Table 1

Results of branch‐site model for SARS‐CoV‐2

Model	Ln L	Parameter estimation					model comparison	LRT p values	Positive selection site
Model A	−157,747.79203	Site type	0	1	2a	2b	Model A vs. Model A null	0.00000	498V* 1,039M* 1,125V* 1,183 E* 1,388N* 1,592R* 1,968P* 2,001S* 2,020S* 2,169S* 2,233S* 2,257S* 2,360S*
		Site Ratio f	0.41494	0.35320	0.12525	0.10661
		Background branch ω0	0.31558	1.00000	0.31558	1.00000
		Detection of branch ω1	0.31558	1.00000	65.38486	65.38486
Model A null	−157,831.26325	1

RDP4 and SimPlot were used to detect the recombination event of the 15 coronaviruses. However, the results showed that no statistically significant recombination event was found.

Harmful mutation detection in S protein amino acids

We found that there were 101 nonsynonymous mutations in all the 69 098 protein sequences. The results showed that there were seven harmful mutations (Table 2), while the remaining 94 were neutral.

Table 2

Seven harmful mutations detected in the S protein

Mutation	Score	Prediction (critical = −2.5)
P589S	−3.966	Deleterious
T716I	−3.293	Deleterious
D936Y	−2.602	Deleterious
S939F	−3.094	Deleterious
P1162S	−2.722	Deleterious
C1236F	−4.061	Deleterious
C1250F	−5.057	Deleterious

DISCUSSION

SARS‐CoV‐2 has evolved into thousands of variants worldwide. Although the virus might mutate a lot, only a few could cause serious harm to humans. Compared with other coronaviruses and RNA viruses, the mutation rate of SARS‐CoV‐2 (6.677 × 10−4) is at a medium level. Moreover, the mutation rate gradually decreased during the subsequent infection process. , , The reason for this may be that SARS‐CoV‐2 has its own copy “proofreading mechanism”, which can correct some errors that may occur during the copying process, leading to the decrease of the mutation rate of SARS‐CoV‐2. On the other hand, harmful mutations may lead to a certain degree of protein structure and functional changes, ultimately affecting the reproduction of SARS‐CoV‐2. Under the pressure of natural selection and with the accumulation of harmful mutations, the number of viruses will gradually decrease and even mutational meltdown may occur, leading to population extinction, , with a result of the relatively low evolution rate of SARS‐CoV‐2 (Table 3).

Table 3

The nucleotide mutation rate (substitutions per site per year) of different RNA virus

Group	Family	Virus	Mutation rate	Reference
ss(+)RNA	Coronaviridae	SARS‐CoV‐2	8 × 10−4	19
ss(+)RNA	Coronaviridae	SARS	3.01 × 10−3	17
ss(+)RNA	Coronaviridae	MERS‐CoV	1.12 × 10−3	15
ss(+)RNA	Coronaviridae	HCoV‐OC43	1.06 × 10−4	41
ss(+)RNA	Coronaviridae	HCoV‐229E	3.28 × 10−4	42
ss(+)RNA	Coronaviridae	Avian coronavirus	2.40 × 10−4	43
ss(+)RNA	Coronaviridae	Bovine coronavirus	5.37 × 10−4	44
ss(+)RNA	Filoviridae	EBOV	1.23 × 10−3	45
ss(+)RNA	Picornaviridae	Hepatitis A virus	9.76 × 10−4	46
ss(+)RNA	Flaviviridae	Hepatitis C virus	1.39 × 10−3	47
ss(−)RNA	Orthomyxoviridae	Influenza A virus	3.15 × 10−3	48
ss(+)RNA	Flaviviridae	Dengue virus	6.50 × 10−4	49
ss(+)RNA	Picornaviridae	Human enterovirus A	5.53 × 10−3	50
ss(+)RNA	Picornaviridae	Human enterovirus B	5.27 × 10−3	50
ss(+)RNA	Picornaviridae	Poliovirus 1	1.17 × 10−2	50
ss(−)RNA	Paramyxoviridae	Measles virus	6.02 × 10−4	51
ss(−)RNA	Rhabdoviridae	Rabies virus	3.32 × 10−4	52
dsRNA	Reoviridae	Human rotavirus A	1.87 × 10−3	53

However, given that SARS‐CoV‐2 has not yet been effectively controlled in the United States, we cannot rule out the possibility that more new mutations will appear in the United States. Compared with other coronaviruses, 13 positive selection sites were detected in the SARS‐CoV‐2 genome, indicating that SARS‐CoV‐2 has a special evolution pattern.

In addition, nonsynonymous mutations were found in many samples of the S protein, with a harmful mutation proportion of 6.93%. These harmful mutations may affect the structure and function of the S proteins. Once the amino acid of RBD of S protein is mutated, the binding affinity with the human ACE2 receptor may change, which may cause the increase of the ability to infect humans and also make the existing vaccines ineffective. Global mutation data collected from the GISAID database revealed that mutations occurred at almost every site of the S gene. However, further experiments are needed to verify whether these mutations would indeed affect the function of the S protein.

Viruses are the masters of evolution, which create new variants by mutating and recombining in an unpredictable way during each replication cycle. From the perspective of biological evolution, the existence of viruses is a natural selection pressure for human beings. The human body's immune response will produce a certain degree of adaptability, which will promote the development of human beings in a direction that is more conducive to human survival. However, the evolution of the human immune system and various interventions may also promote the adaptive mutations of the virus. Usually, viruses become “mild” as they circulate in the same host because high pathogenicity will lead to the death of the host, resulting in loss of transmission and reproduction, which is not conducive to the survival and reproduction of the virus itself. As mentioned above, the current mutation rate of SARS‐CoV‐2 in the United States is lower than that at the beginning of the outbreak. Obviously, this view is in line with our expectations.

Experimental results showed that some variants of SARS‐CoV‐2 strains increased the infectivity, such as variants D614G, S477N, and N439K, and these variants had higher transmission ability and faster replication speed than that of the original viruses, however, their pathogenicity did not increase. , , ,

At present, more infectious variants of SARS‐CoV‐2 appearing in the world include the first D614G variant B.1.5‐B.1.72, British variant B.1.1.7, South African variant B.1.351, Brazilian variant P.1, and Philippine variant P.3, as well as the two native varieties of the United States, the California variety B.1.429/B.1.427 and New York variety B.1.526. The spreading power, virulence, and immune evasion ability of new variants have been increasing. Among them, the British variant B.1.1.7 is a virus that has a relatively large impact on the population. In December 2020, the new strain B.1.1.7 was discovered in the UK for the first time, and its S gene sequence accumulated 16 nucleotide mutations, resulting in 10 amino acid site changes (H69del, V70del, Y144del, N501Y, A570D, D614G, P681H, T716I, S982A, and D1118H). According to the data from PANGO lineages, B.1.1.7 mutant strain has been found in more than 90 countries. Studies have shown that the infectiousness of this mutant strain is more than 50% higher than that of the current prevailing strains. As of April 10, the strain has caused 20,915 confirmed infections in the United States and has become the main type of transmission in the United States. The California variant B.1.429/B.1.427 discovered in July 2020 has a mutation site of L452R, which leads to a 20% increase in its transmission power, and exhibits moderate immune evasion leading to a fairly fast transmission rate. The New York variant B.1.526 discovered in November 2020, whose mutation site is E484K or S477N, also exhibits moderate immune evasion.

At present, the main variants concerned in the United States include B1.1.7, B1.351, P.1, B.1.429/B.1.427. There is evidence that these variants could lead to increased infectivity, more severe disease, reduced effectiveness of treatments or vaccines, or diagnostic detection failures. , , At present, the proportion of B1.1.7 has reached 44.1% and shows a significant trend of continuing to increase. Taking New York as an example, B1.1.7 accounted for 11.9% in February, 26.2% in mid‐March, and 28.2% at the end of March. At the same time, the proportion of the New York native variant B.1.526 (E484K) is still as high as 27.9%. These different variants have greatly affected the effectiveness of the vaccine.

It is worth noting that the B.1.617 variant strain first appeared in India and has now become one of the most popular strains in the United States, and B.1.617 has a 110% higher affinity for human ACE2 receptor, making it highly infectious. ,

In conclusion, we clarify the evolutionary characteristics of SARS‐CoV‐2 in the United States, providing a scientific basis for future surveillance and prevention of virus variants. To control SARS‐CoV‐2 and to restore people's normal life activities as soon as possible, it is necessary to continue to monitor specific mutations, which is still of great significance for further in‐depth study of SARS‐CoV‐2 and the evaluation of the effectiveness of existing vaccines.

CONFLICT OF INTERESTS

The authors declare that there are no conflicts of interest.

AUTHOR CONTRIBUTIONS

Sihua Peng conceived and designed the study. Shihang Wang, Jingying Zhao, Yin Zheng, and Xiaoyu Liu collected data and prepared the datasets. Shihang Wang, Xuanyu Xu, and Cai Wei participated in phylogenetic analyses. Shihang Wang and Xuanyu Xu drafted the manuscript. Xiaomin Zeng, Wenliang Yuan, and Sihua Peng provided critical comments. All authors contributed to manuscript revision, read, and approved the submitted version.

Supporting information.

Click here for additional data file.