The ongoing outbreak of coronavirus disease 2019 (COVID-19), caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), has led to tremendous human mortality and economic hardship. As of 31 July 2020, over 17,106,007 confirmed COVID-19 cases have been reported worldwide and 668,910 deaths have occurred from the disease [1
]. To mitigate this devastating pandemic, we have to control its spread by sufficient testing, social distancing, contact tracking, and developing effective diagnosis tools, efficacious antiviral drugs, antibody therapies, and preventive vaccines.
SARS-CoV-2 is a positive-sense, single-strand RNA virus that belongs to the beta coronavirus genus [2
]. It has a genome size of 29.82 kb, which encodes multiple non-structural and structural proteins. The ORF1ab encode non-structural proteins for RNA replication and transcription. The downstream regions of the genome encode structural proteins, including the spike (S) protein, the nucleocapsid (N) protein, the envelope (E) protein, and the membrane (M) protein. All of the four major structural proteins are required to produce a structurally complete viral particle. The S protein mediates viral attachment to the host angiotensin-converting enzyme 2 (ACE2) receptor and subsequent fusion between the viral and host cell membranes aided by transmembrane serine protease 2 (TMPRSS2) to allow the entry of viruses into the host cell [3
]. The nucleocapsid (N) protein, one of the most abundant viral proteins, binds to the RNA genome and is involved in replication processes, assembly, and host cellular response during viral infection [6
Mutagenesis is a basic biological process that changes the genetic information of organisms. As a primary source of many kinds of cancers and heritable diseases, mutagenesis may be fearful but is a driving force for natural evolution [7
]. Although viruses are not organisms per se, they are at the edge of life. Our real-time interactive SARS-CoV-2 Mutation Tracker (https://users.math.msu.edu/users/weig/SARS-CoV-2_Mutation_Tracker.html
) shows that over 15,000 mutations have occurred on SARS-CoV-2 [9
]. More than 1700 mutations on the S protein gene have a significant impact on SARS-CoV-2 infectivity [10
]. These mutations should be put into the perspective that COVID-19 has globally spread. The geographical and demographical diversity of the viral transmission and exogenous and endogenous genotoxins exposure have stimulated SARS-CoV-2 mutations. If we consider the average number of mutations per genome, SARS-CoV-2 is mutating slower than other viruses, such as the flu and common cold viruses [13
]. This is because SARS-CoV-2 belongs to the coronaviridae family and the Nidovirales order, which has a genetic proofreading mechanism in its replication achieved by an enzyme called non-structure protein 14 (NSP14) in synergy with NSP12, i.e., RNA-dependent RNA polymerase (RdRp) [15
]. As a result, SARS-CoV-2 has a relatively high fidelity in its transcription and replication processes. In general, Coronavirus mutations are created from three major sources, namely, random errors in replication, such as genetic drift and spontaneous genotoxins, viral replication proofreading and defective repair mechanisms, and host immune responses, such as gene editing [11
]. Genotyping tracks mutations, overpopulation, space, and time, while also providing a method to understand the molecular mechanism of SARS-CoV-2 proteins, protein–protein interactions, and their synergy with host cell proteins, enzymes, and signaling pathways.
The studies of SARS-CoV genomes have, to date, predominantly focused on understanding genome mutation variants, implications in virus transmissions [18
], and ramifications on the development of diagnostics [9
], vaccines [21
], antibodies [22
], and drugs [21
]. Although it is difficult to determine the detailed mechanism of every specific mutation, early work on a few initial SARS-CoV-2 strains in Wuhan, China, revealed that hypermutations C > T are most likely resulted from the APOBEC (apolipoprotein B mRNA editing enzyme, catalytic polypeptide-like) deamination in RNA editing [23
]. Simmonds noted that a large proportion of C > T mutations in a host APOBEC-like context in the initial months of the pandemic provide evidence for a potent host-driven antiviral editing mechanism against SARS-CoV-2 [24
]. An accumulation of C > T mutations in SARS-CoV-2 variants was also reported in [25
]. In the standard genetic code, all three stop codons, TAA, TAG, and TGA, involve T but not C. Therefore, the gene-editing imposed C > T mutations will have a high possibility of terminating the translation of viral proteins, which undermines viral functions and survivability. To be noted, as stated in [25
], the C > T mutation can increase the frequency of codons for hydrophobic amino acids, which may have an effect on the properties of SARS-CoV-2 proteins. The spontaneous C > T transitions and APOBEC deamination can lead to cancer in humans. There are two well-known deaminase RNA editing mechanisms in human cells: the APOBEC [26
] and the adenosine deaminases acting on RNA (ADAR) [27
]. The APOBEC enzymes deaminate cytosines into uracils (C > U) on single-stranded nucleic acids (ssDNA or ssRNA). It is well established that the human genome encodes activation-induced cytidine deaminases (AIDs) and several homologous APOBEC cytidine deaminases that function in innate immunity as well as in RNA editing [28
]. In both innate and adaptive immunity, AID and APOBEC cytidine deaminases modulate immune responses by mutating specific nucleic acid sequences of hosts and pathogens. The ADAR enzymes deaminate adenines into inosines (A-to-I) and result in A > G mutation. The significance of A-to-I editing is its abundance in both host and viral RNAs. ADAR enzymes play important roles during viral infections. They can have either a proviral or an antiviral consequence, depending on the virus–host combination [30
The APOBEC family proteins play critical functional roles within the adaptive and innate immune system, which is involved soon after the infection [32
]. Therefore, the higher ratio of C > T mutations may indicate the strong capacity of the host immune system. However, a strong immune response is a double-edged sword. On the one hand, it may help host cells to defeat the virus more efficiently. On the other hand, it can result in a “cytokine storm”, which is a key cause of the death of COVID-19 patients by the exponential growth of inflammation and organ damage [33
In this work, we analyze a large volume of single-nucleotide polymorphisms (SNPs) found in 33,693 complete SARS-CoV-2 genome isolates globally. By analyzing the distribution of 12 SNP types, we notice that the ratio of C > T mutations is predominately higher than that of the other types of mutations, indicating that hypermutation C > T may result from extensive host RNA editing, i.e., the APOBEC deamination. Additionally, we investigate the distribution of 12 SNP types in different age groups, gender groups, and geographic locations to understand whether these hypermutations have a age/gender/demographic preference. Moreover, we provide deep insights into the mutation motif and hot-spot patterns from 13,833 single mutations decoded from 33,693 complete SARS-CoV-2 genome sequences, revealing mutational signatures and preferred genetic environments. Finally, we hypothesize that virus genomes evolve through host innate immune response imposed gene editing, i.e., C > T, and virus protective mechanism-installed defective revisionary mutations, T > C. As a result, both C > T and T > C mutation ratios are usually high. We show that the ratio of C > T to T > C mutations is higher than the unity in the forward viral evolution, which suggests the master and slave relationship between host gene editing and virus protective mechanism. Therefore, we propose that the C > T to T > C ratio being higher than the unity (>1) is an indication of the forward viral evolution direction.
To reveal that C > T and A > G mutations are driven by RNA-APOBEC and RNA-ADAR editing, we first analyzed 33,693 complete SARS-CoV-2 genome sequences, and a total of 13,833 single mutations were found as of 31 July 2020. To be noted, 13,833 single mutations were unique mutations, i.e., the same mutation appearing in different SARS-CoV-2 isolates is only counted once. If we count the same mutation in different SARS-CoV-2 isolates repeatedly according to their frequency, then all of the mutations that are detected in the 33,693 complete SARS-CoV-2 genome sequences are called non-unique mutations. With the reference sequence of SARS-CoV-2 genome collected on 5 January 2020 [2
], we calculated the proportion of 12 SNP types (i.e., A > T, A > C, A > G, T > A, T > C, T > G, C > T, C > A, C > G, G > T, G > C, G > A) worldwide. The unusually high ratios of C > T and A > G mutations indicate that RNA-APOBEC editing and RNA-ADAR editing are involved in the host immune response to SARS-CoV-2 infection [28
]. Additionally, to understand gene-editing preference, we investigated the distribution of 12 SNP types of mutations in different countries/regions, age groups, and gender groups. Furthermore, we decoded mutation motifs from the 2-mer and 3-mer sequence contexts to survey the hot-spot patterns and mutational signatures driven by gene-editing. Moreover, we analyzed the proportion of 12 SNP types among SARS-CoV, Bat-SL-BM48-31, Bat-SL-CoVZC45, Bat-SL-RaTG13, and SARS-CoV-2. We discovered that the viral evolution order can be determined by the ratios of C > T/T > C. These results are presented in following subsections.
3.1. Host Immune Response to SARS-CoV-2 Infection with Gene Editing
3.1.1. Global Analysis
illustrates the proportion of 12 SNP types of SARS-CoV-2 (i.e., A > T, A > C, A > G, T > A, T > C, T > G, C > T, C > A, C > G, G > T, G > C, G > A) in the global. Here we only consider the unique SNPs.
First, it can be seen that not all SARS-CoV-2 mutations are created equal. Mutation C > G only accounts for 1.25%. A few other mutation types, G > C, T > G, T > A, and A > C, are not frequent either. If mutations are random, each mutation should have a ratio of 8.3% on average. It can be seen that C > T owns the largest proportion (24.06%), which is much higher than the average ratio. Applying the Z-test, the p-value is , which means our deduction is of statistical significance. Therefore, the hypermutation C > T must be driven by additional mechanisms. It is known that host RNA-APOBEC editing leads to excessive C > T transitions.
Moreover, the second most frequent mutation type is A > G transition. Its A > G ratio of 14.87% is much higher than the average ratio of 8.3%, and the p-value from the Z-test implies that our deduction passes the hypothesis test. Therefore, RNA-ADAR editing is also involved in the host immune response. Although the high ratios of C > T and A > G reveal that the immune system is combating SARS-CoV-2 by two deaminase RNA editing mechanisms, the relatively high ratios of the reversed mutations T > C and G > A also indicate that SARS-CoV-2 fights the gene editing using its defective proofreading and repairing mechanisms.
Finally, it is well-known that mutations can be classified into four transition types (i.e., A > G, G > A, C > T, and T > C) and eight transversion types. Table 1
shows that all transition types have relatively high ratios, whereas all transversion types, except for G > T, have relatively low ratios. This is due to the fact that it is easier to substitute a single-ring nucleotide structure for another single-ring nucleotide structure than to substitute a double-ring nucleotide for a single-ring nucleotide. Additionally, transitions are more likely to result in silent mutations. Therefore, transversions can be more destructive to viral genomes.
3.1.2. Age Analysis
illustrates the distribution of 12 SNP types among unique SNPs in SASR-CoV-2 genome isolates from different age groups. In general, with the increase in age, the ratio of C > T gradually increased. Here, 42.1% C > T mutations were detected in patients who are older than 90 years old, indicating that gene-editing is more active in elderly patients than young patients. We applied the Z-test to the data of patients older than 90 years old and patients in other age groups, and all the resulting p
-values were less than
. Therefore, our hypothesis is statistically significant. However, the severe COVID-19 cases may be due to the immune systems’ heightened response. When SARS-CoV-2 infects a host cell, a set of proteins called cytokines will be released from a broad range of cells (mainly immune cells). Cytokines are involved in the immune response to produce more immune cells and recruit them to the sites of inflammation in order to fight against the viral infection. In turn, more cytokines can be released from the immune cells. This positive feedback loop will result in a “cytokine storm”, which can beget the exponential growth of inflammation, trigger apoptosis, and lead to organ damage [33
]. Therefore, we hypothesize that if the immune system overreacts to the invading pathogens, it is more likely to cause the cytokine storm and aggravate the condition of COVID-19 patients. It can be seen in Figure 1
, where patients who are older than 80 years old have more C > T mutations compared to other age groups. This result reveals that APOBEC3 gene editing is more active with older people. Consequently, the cytokine storm may happen more frequently in older people than it does in younger people. This might be one of the main causes of the high COVID-19 fatality for the elderly. Age-related mutagenesis, i.e., C > T transition, is known to cause more cancer diagnostics in the elderly [39
Notably, the SARS-CoV-2 samples from children under five years old have a relatively high ratio of C > T mutations (39.6%). The p-value from the Z-test is less than , which also implies the statistical significance of a relatively high ratio of C > T in patients under five years old. Moreover, the reversed mutation type T > C for samples from children under five years old and adults older than 90 years old has the third-largest ratio. In other age groups, T > C has the fourth-largest ratio. As demonstrated before, the reversed mutation T > C may reveal that SARS-CoV-2 is capable of fighting back against the host immune system. Therefore, we deduce that SARS-CoV-2 will fiercely counter-attack the immune system in children under five and adults older than 90.
Our result reveals that the immune systems of children under five years old are less well-developed and weaker than those of adults. They have to fight harder when SARS-CoV-2 infects. This result suggests that children under five are at risk of COVID-19. However, the long-term health consequence of young children’s unusual response to SARS-CoV-2 infection remains to be further studied.
3.1.3. Gender Analysis
shows the distribution of 12 SNP types in SARS-CoV-2 genome isolates globally from two gender groups. The ratio of C > T mutations in females is slightly higher in males. We also apply a simple Z-test, and the p
, which is smaller than the significance level
. Therefore, the slightly higher C > T ratios in females matches the finding that women have a stronger immune response than men [40
]. Moreover, Figure 3
and Figure 4
depict the distribution of 12 SNP types in different age groups among female and male patients. Overall, the proportion of C > T mutations in the SARS-CoV-2 genomes from females is higher than the C > T proportion in the SARS-CoV-2 genomes from males, except for ages between 6 and 19 and ages older than 90. Therefore, we can deduce that the RNA editing has an age and gender preference, and is therefore more likely to happen or become stronger for females who are older than 90 years old or under 5 years old.
3.1.4. Geographic Analysis
In this section, we analyze the distribution of SARS-CoV-2 mutations in different countries and regions. Limited by the number of complete genome sequences submitted to GISAID that have appropriate labels, we only analyzed the countries with more than 1000 labeled sequences to maintain statistical significance. Table 2
lists the total number of SARS-CoV-2 sequences in the United Kingdom, United States, Australia, and India. The number of sequences with age and gender information is given in Table 2
illustrates the C > T ratios in the SARS-CoV-2 genome isolates from different age groups in the United Kingdom, United States, Australia, and India, respectively. We can see that the SARS-CoV-2 genome isolates from the United Kingdom patients have the highest ratio of C > T compared to those from the other three countries. It is interesting to note that the SARS-CoV-2 genome isolates from the patients older than 80 years old from the United Kingdom and Australia have less C > T mutations, which is not consistent with the global pattern. The distribution of 12 SNP types in the SARS-CoV-2 genome isolates from different age groups in the United Kingdom, United States, Australia, and India can be found in the Supporting Information
illustrates the distribution of 12 SNP types in six continents. The SARS-CoV-2 genome isolates from Europe, Asia, and North America patients have a relatively low C > T mutation ratio (less than 35%), while the reversed T > C mutation ratio is relatively high (greater than 10%). When employing the Z-test between the data in the world and in Europe, Asia, and North America, the p
-values are all less than
, which means that the C > T and T > C ratios in Europe, Asia, and North America are not at the same level as the global ratio. On the contrary, South America, Oceania, and Africa have higher C > T ratios but a lower T > C ratio. It worth noting that the C > T mutation ratios in the SARS-CoV-2 genome isolates from Oceania and Africa are more than 10% higher than those of Asia, Europe, and North America. This result indicates that the APOBEC editing may be more active, and the counterattack of SARS-CoV-2 might be weakened by the strong immune response in the populations of Oceania and Africa.
African Americans, as an ethnic group of Americans with total or partial ancestry from Africa, are genetically associated with Africans. There have been many concerns about the fact that they are disproportionately affected by COVID-19 (https://www.cdc.gov/coronavirus/2019-ncov/community/health-equity/race-ethnicity.html
). The present finding indicates that the immune systems of African Americans may also overreact to SARS-CoV-2 infection by excessive gene editing.
Another interesting issue is that A > G mutation ratios of genome isolates from Oceania and Africa are very low (<9%). In contrast, the A > G mutation ratios of genome isolates from other regions are significantly higher (>11%). The Z-test was also employed to prove that our deduction is statistically significant. For example, when we extracted the number of A > G mutations in Oceania and Asia for calculation, the p-value was , which is less than the significance level . Similarly, we repeated the procedure. For patients in Oceania and North America, Oceania and South America, Oceania and Europe, all the p-values are less than the significance level . These results indicate that Asia and Europe populations may have adopted significantly different genetic and molecular mechanisms in their immune response to viral infection compared to those of Oceania and Africa. Further studies are required to fully understand these differences.
3.2. The SNP Preferences on Sequence Contexts
The mutation preferences in sequence contexts may be used to predict the mutational signatures from genome sequences. Despite numerous studies of the mutation contexts in APOBEC editing inhuman cells, little is known of the mutation contexts in the SARS-CoV-2 genome. As we have a large number of SNP mutations from SARS-CoV-2 genomes, here, we discuss the mutation frequencies from 2-mer and 3-mer sequence contexts. We present 4-mer sequence contexts in the Supporting Information
In general, the patterns discussed in this section are consistent with those presented in Section 3.1.1
. However, this section offers more detailed information about mutational signatures.
For mutation motifs of SNPs at the first position of 2-mers Figure 6
a, we observe that motif 2-mer C
W (where W is either A or T) for the C > T mutation is the predominant context. Similarly, for mutation motifs of the SNPs at the second position of 2-mers Figure 6
b, motif 2-mer WC
for C > T mutation is the predominant context. These results are consistent with the previous study that TC
W contexts (where W = A or T) are predominantly caused by the APOBEC-catalyzed deamination of cytosine (C) to thymine (T) or uracil (U) in human cancer cells [42
For the SNPs at the first position of 3-mers (A
NN or T
NN) (Figure 7
a), we observe the following mutation patterns.
ANN has high A > G mutations;
TNN has a high frequency in T > C mutations.
For the SNPs at the first position of 3-mers (C
NN or G
NN) shown in Figure 7
b, we observe the following mutation patterns.
CNN has a high frequency in C > T mutations;
GGA and GGT has a high frequency in G > A mutations;
GAA has a relatively high frequency in G > A mutations;
GGW (where W is A or T) has relatively high frequency in G > T mutations.
For the SNPs at the second position of the 3-mers (NA
N or NT
N), as shown in Figure 8
a, we observe the following mutation patterns.
NAN has a high frequency in A > G mutations;
NTN has a high frequency in T > C mutations;
The A > C mutation also has a larger proportion in AAW (where W is A or T).
For the SNPs at the second position of the 3-mers (NC
N or NG
N), as shown in Figure 8
b, we observe the following mutation patterns.
WGN (where W is A or T) has a G > T dominated mutation except for AGG;
SGN (where S is G or C) has G > A dominated mutations;
AGG has high G > A mutations;
Characteristic combinations SCG (where S is G or C) are stable and only a few C > G mutations are detected;
Characteristic combinations GGS (where S is G or C) are stable, only a few G > C mutations are detected.
For the SNPs at the third position of 3-mers (NNA
) a shown in Figure 9
a, we observe the following mutation patterns.
A > G mutation has a high frequency in NNA;
T > C mutation has a high frequency in NNT;
T > C mutation is dominated in NGT and only a few T > A and T > G are found in the sequence context of NGT.
For the SNPs at the third position of 3-mers (NNC
) as shown in Figure 9
b, we observe the following mutation patterns.
NNC has a high frequency in C > T mutations;
G > T mutation has a high frequency in NNG;
G > A is also highly expressed in the sequence context of NNG;
Characteristic combinations CGC are stable and the mutations on these patterns are most likely to be C > T transitions.
3.3. Coronavirus Evolution
It is reasonable to assume that the five coronaviruses, SARS-CoV (2003) [43
], Bat-SL-BM48-31 (2008) [44
], Bat-SL-CoVZC45 (2017) [45
], Bat-SL-RaTG13 (2013) [46
] and SARS-CoV-2 (2019) [2
], are of the same origin but differ from each other by their evolutionary stages. The data collection date of Bat-SL-RaTG13 (2013) was denoted as 24 July 2013, while the data were not uploaded to the GISIAD database until 27 January 2020. Figure 10
shows the mutation ratio among these five genomes. First, similar to the SARS-CoV-2 mutations listed in Table 1
, four transition types (i.e., A > C, C > A, C > T, and T > C) still have high mutation ratios. Notably, the C > T type has the highest ratio, indicating that host immune response still plays the major role. However, transversion type G > T is not as important as in the SARS-CoV-2 mutations discussed earlier. Nonetheless, transversion types A > T and T > A appear on the top six mutation types.
We hypothesize that gene editing via APOBEC (C > T) and ADAR (A > G) is a driving force for RNA viral evolution, as shown in Table 1
. Viruses may fight the host immune response with either defective repair or reversed mutations (T > C) within survived isolates. Therefore, the T > C mutation rate would decrease during evolution. We are interested in not only the C > T transition ratio, but also the ratio of C > T over T > C, the reversed transitions. From Figure 10
, we can calculate the C > T over T > C ratios among different SARS-like species which are listed in Table 4
It is seen that viral evolution order may be determined by the T > C over T > C ratio. Through this analysis, we have the following evolution order for the aforementioned coronaviruses: SARS-CoV (2003) → Bat-SL-BM48-31 (2008) → Bat-SL-CoVZC45 (2017) → Bat-SL-RaTG13 (2013) → SARS-CoV-2 (2019) → 33693 SARS-CoV-2 genome isolates (2020). Here, we have one reversed order between Bat-SL-CoVZC45 (2017) → Bat-SL-RaTG13 (2013). This may happen for a few reasons. First, these coronaviruses may not be of the same origin. Second, the data collection date may not be accurate. The sequence of Bat-SL-RaTG13 (2013) was not uploaded until 2020. Finally, our method may admit a few counterexamples.