key: cord-0899606-0swnxe2a
authors: Huang, Wanyi; Guo, Yaqiong; Li, Na; Feng, Yaoyu; Xiao, Lihua
title: Codon usage analysis of zoonotic coronaviruses reveals lower adaptation to humans by SARS-CoV-2
date: 2021-01-28
journal: Infect Genet Evol
DOI: 10.1016/j.meegid.2021.104736
sha: 90ec5fbf9edd270f4357ddb9d08bbbe110138b93
doc_id: 899606
cord_uid: 0swnxe2a

Since 2002, the world has witnessed major outbreaks of acute respiratory illness by three zoonotic coronaviruses (CoVs), which differ from each other in pathogenicity. Reasons for the lower pathogenicity of SARS-CoV-2 than the other two zoonotic coronaviruses, SARS-CoV and MERS-CoV, are not well understood. We herein compared the codon usage patterns of the three zoonotic CoVs causing severe acute respiratory syndromes and four human-specific CoVs (NL63, 229E, OC43, and HKU1) causing mild diseases. We found that the seven viruses have different codon usages, with SARS-CoV-2 having the lowest effective number of codons (ENC) among the zoonotic CoVs. Human codon adaptation index (CAI) analysis revealed that the CAI value of SARS-CoV-2 is the lowest among the zoonotic CoVs. The ENC and CAI values of SARS-CoV-2 were more similar to those of the less-pathogenic human-specific CoVs. To further investigate adaptive evolution within SARS-CoV-2, we examined codon usage patterns in 3573 genomes of SARS-CoV-2 collected over the initial 4 months of the pandemic. We showed that the ENC values and the CAI values of SARS-CoV-2 were decreasing over the period. The low ENC and CAI values could be responsible for the lower pathogenicity of SARS-CoV-2. While mutational pressure appears to shape codon adaptation in the overall genomes of SARS-CoV-2 and other zoonotic CoVs, the E gene of SARS-CoV-2, which has the highest codon usage bias, appears to be under strong natural selection. Data from the study contribute to our understanding of the pathogenicity and evolution of SARS-CoV-2 in humans.

maintaining the amino acid sequences, leading to changes in viral replication (Chen et al., 2020) . Effective number of codons (ENC) and codon adaptation index (CAI) analyses are frequently used to evaluate the codon usage patterns of viruses (Baha et al., 2019; Belalov and Lukashev, 2013; Khandia et al., 2019; G. Li et al., 2018; Stoletzki and Eyre-Walker, 2007; Zhang et al., 2019) . Between them, ENC values are measurements of codon preference (Comeron and Aguadé, 1998) , while CAI is an indirect measurement of the likely expression efficiency of viral genes in the host (Carbone et al., 2003; Chen et al., 2020) .

In this study, we performed comparative analyses of codon usage patterns among the three zoonotic CoVs and the four human-specific CoVs to understand the differences in pathogenicity. We also examined the evolution of codon usage in SARS-CoV-2 genomes collected over a 4-month period during the early pandemic.

A total of 122 complete CoV genomes (15 from SARS-CoV, 15 from SARS-CoV-2, 25 from MERS-CoV, 15 from HCoV-OC43, 15 from HCoV-HKU1, 15 from HCoV-229E, 15 from HCoV-NL63, 15 from Bat-CoV, and 2 from SARS-CoV-2-related CoV) were downloaded from the National Center for Biotechnology Information (https://www.ncbi.nlm.nih.gov/genbank/). After removing low-quality sequences, 86 CoV genomes were used in analyses.

A total of 17,440 SARS-CoV-2 genomes flagged as "collection date from December 1 2019 to April 30 2020" were downloaded from the EpiCoV database of the GISAID Initiative (https://www.epicov.org) as of May 10 2020. After filtering to exclude those with lengths less than 29,000 bp, those with Ns or degenerate bases, and those derived from animals, a final dataset of 3,573 genomes was used in downstream analyses.

Reference genomes of human and zoonotic CoVs were obtained from GenBank (https://www.ncbi.nlm.nih.gov/genbank/). The open reading frames of the genomes were predicted using Geneious v11.1.5 (https://www.geneious.com/) and annotated using Blastn v2.10.1 (https://blast.ncbi.nlm.nih.gov/Blast.cgi).

The CoV genomes were aligned using MUSCLE v3.8.31 (https://www.ebi.ac.uk/). The substitution model was selected using jModelTest v2.1.10 (Posada, 2008) , based on values from the Akaike Information Criterion (AIC), the Bayesian Information Criterion (BIC), and hierarchical likelihood-ratio tests (hLRTs). A maximum likelihood (ML) tree was reconstructed using PhyML (Guindon et al., 2010) , with the GTR + GAMMA I + I model and bootstrap value of 1,000.

The frequency of each nucleotide (A, U, G, C), AU, and GC in the genes or J o u r n a l P r e -p r o o f Journal Pre-proof genomes was calculated using BioEdit v7.1.3.0 (https://bioedit.software.informer.com/). The nucleotide frequency of the third position of synonymous codons (A 3 %, U 3 %, G 3 %, C 3 %) was calculated using CodonW (http://codonw.sourceforge.net/). The GC content at the first (GC 1S ), second (GC 2S ), and third codon positions (GC 3S ) was calculated using Emboss explorer (http://codonw.sourceforge.net/), with the GC 12S as the mean of GC 1S and GC 2S . The three termination codons (UAA, UGA, UAG) and the codons AUG and UGG were excluded from these analyses.

To assess the codon usage patterns without the effect of sequence length, RSCU values for 59 codons (excluding UAA, UGA, UAG, AUG, and UGG) were calculated using DAMBE v7.2.43 (http://dambe.bio.uottawa.ca/). An RSCU value of 1 indicates that the codon is used equally, while codons with RSCU values of >1.6 and <0.6 are considered as over-represented and under-represented, respectively (Sharp and Li, 1986) . The principal component analysis (PCA) implemented in TBtools v0.66831 (http://www.tbtools.com/) was used in the analysis of the RSCU data based on a 59-dimension vector.

The ENC was calculated using the CodonW software v1.4.4 (http://codonw.sourceforge.net/). The ENC values range from 20 (only one synonymous codon is used per amino acid, showing an extreme codon usage bias) to 61 (all synonymous codons are equally used showing no bias) (Wong et al., 2010) . In general, ENC values of less than 35 indicate strong codon usage bias (Comeron and Aguadé, 1998) . The ENC values generated were plotted against GC 3S . If the codon usage bias is merely constrained by mutational pressure, the ENC values would lie on or around the expected standard curve. Values below the expected curve, on the other hand, indicate that natural selection could play a role in shaping the codon usage (Wong et al., 2010) . Plots of the expected ENC values against GC 3S values were generated as described (Zhang et al., 2019) . eight clades of CoVs analyzed in the study (Figure 1A) , Clade 3 (HCoV-HKU1) had the lowest GC content (32.6%) and the lowest GC content at the third codon position (GC 3 ) (18.9%). In contrast, the highest values of both indices (GC = 41.7% and GC 3 = 36.3%) were seen in Clade 7 (Bat-CoV). Specifically among the three zoonotic CoVs, SARS-CoV-2 (Clade 6) had the lowest GC content (39.2%) and the lowest GC 3 content (30.1%) ( Table 1) , while MERS-CoV (Clade 5) had the highest values of both indices (GC = 41.6% and GC 3 = 35.6%). Moreover, GC contents of other clades were 35.9%, 37.1%, 39.1%, and 40.9% for Clade 2 (HCoV-NL63), Clade 4 (HCoV-OC43), Clade 1 (HCoV-229E), and Clade 8 (SARS-CoV), respectively. In comparison, the GC 3 contents of those clades were 24.7%, 28.3%, 32.4%, and 33.8%, respectively.

In the ML analysis, the CoVs under analysis formed eight clades with high bootstrap support: Clade 1 (HCoV-229E), Clade 2 (HCoV-NL63), Clade 3 (HCoV-HKU1), Clade 4 (HCoV-OC43), Clade 5 (MERS-CoV), Clade 6 (SARS-CoV-2), Clade 7 (Bat-CoV), and Clade 8 (SARS-CoV) ( Figure 1A ). In addition, the ML tree showed that the genome KT253327 was genetically related to HCoV-229E, the genome KY417150 was genetically related to SARS-CoV, and the pangolin CoV EPIISL410721 and the bat CoV MN996532 were genetically related to SARS-CoV-2. Therefore, KT253327 was used as the control for HCoV-229E, KY417150 as the control for SARS-CoV, and EPIISL410721 and KY417150 as controls for SARS-CoV-2 in other analyses in the study. In PCA analysis of the RSCU data (Table S1 ), all clades formed separate clusters together with their controls ( Figure 1B) .

Comparison of the predicted coding regions of the three zoonotic CoVs and four human-specific CoVs showed that they possessed a similar genomic organization ( Figure 1C ). At least six common coding regions were predicted, including 1ab, 1a, S, E, M, and N ( Figure 1C) . The lengths and the order of these genes were similar among the seven CoVs.

In the ENC analysis of the CoV genomes, the values generated were all higher than 35, indicating low codon usage bias by CoVs. Among the eight major phylogenetic clades, the highest ENC value was 52.02 ± 0.30 from the Clade 7 (Bat-CoVs), indicating that these CoVs had the lowest codon usage bias among the lineages. In contrast, the lowest ENC value of human-specific CoVs was 36.27 ± 0.22 from Clade3 (HCoV-HKU1). Among the three zoonotic CoVs, the lowest ENC value was 47.60 ± 0.02 from Clade 6 (SARS-CoV-2), suggesting that SARS-CoV-2 uses a narrower set of synonymous codons (Figure 2A) . This was also the case in ENC analysis of individual genes, especially those from the E gene. In the latter, the ENC value was 42.00 ± 0.00 in Clade 6 (SARS-CoV-2), compared with 61.00 ± 0.00 and 55.71 ± 0.81 in Clades 8 (SARS-CoV) and 5 (MERS-CoV), respectively (Figure 2A ).

To explore whether mutational pressure, natural selection, or they both shaped the codon usage in CoVs, the ENC values were plotted against the GC 3S . The plots J o u r n a l P r e -p r o o f generated with data from full genomes and individual genes showed that the points clustered mostly on or near the expected curves, indicating that mutational pressure was largely responsible for the codon usage bias in CoVs ( Figure 2B ). For the E gene, data from Clade 6 (SARS-CoV-2) and its relative controls (Controls 2 & 3), however, clustered under the expected curve ( Figure 2B) , indicating that natural selection shaped codon usage of this gene in SARS-CoV-2.

We further compared the extent of adaptation to human codon usage by CoVs using CAI analysis. The results obtained suggested that among the eight phylogenetic groups, Clade 5 (MERS-CoV) had the highest CAI value (0.698), followed by Clade 7 (bat-CoVs; 0.690), Clade 8 (SARS-CoV; 0.689), Clade 1 (HCoV-229E; 0.683), Clade 4 (HCoV-OC43; 0.676), Clade 6 (SARS-CoV-2; 0.674), Clade 2 (HCoV-NL63; 0.658), and Clade 3 (HCoV-HKU1; 0.655) ( Figure 2C) . The difference among clades was significant (p < 0.01). Among the three zoonotic CoVs, the lower CAI value in SARS-CoV-2 supports the suggestion that the gene expression of SARS-CoV-2 could be less efficient than that of SARS-CoV and MERS-CoV.

We calculated the ENC values of SARS-CoV-2 isolates collected over a four-month period to evaluate if codon usage changed in SARS-CoV-2. The ENC values decreased gradually during the period, being 47.601 ± 0.014, 47.600 ± 0.015, 47.591 ± 0.016, and 47.588 ± 0.020 for January, February, March, and April 2020, respectively, indicating the codon usage bias of SARS-CoV-2 was increasing slowly over time ( Figure 3A) .

In the longitudinal evaluations of changes in adaptation to human codon usage by SARS-CoV-2, the CAI values of the viral genomes decreased gradually over the four months (p > 0.05), indicating that the likely efficiency of gene expression of SARS-CoV-2 in the human host could be decreasing. This was also the case for most individual viral genes. For example, the CAI values decreased from 0.67200 to 0.67194, 0.67405 to 0.67398, and 0.70507 to 0.70493 for the RdRp, S, and N, genes over the 4-month period, respectively (p < 0.01) ( Figure 3B ). For the M gene, the CAI values also decreased from 0.65106 to 0.65102 over the 4 months, although there was a spike in March when the CAI value reached 0.65115 (p < 0.01). In contrast, CAI values of the E gene remained stable during the first three months and increased from 0.58300 to 0.58305 in April (p = 0.1).

We further evaluated the role of mutational pressure in the adaptive evolution of SARS-CoV-2. The ENC-GC 3S plots generated with data from the entire genome and individual genes showed that the points mostly clustered on or near the expected curves ( Figure S1 ). For the E gene, however, the points were dispersed under the expected curve (Figure S1 ), indicating that unlike other genes, the E gene was further enduring natural selection during the 4-month period.

SARS-CoV-2 is the third zoonotic CoV responsible for a major epidemic after SARS-CoV and MERS-CoV but has spread at unprecedented speed across the globe.

There are four other more common coronaviruses causing mild, self-limiting respiratory symptoms in humans, including HCoV-229E, HCoV-NL63, HCoV-OC43, and HCoV-HKU1. While these CoVs differ from each other in pathogenicity and transmissibility, the reasons for the differences are poorly understood (Chen, 2020) . In this study, we have performed comparative analyses of codon usage among them. We found that the seven viruses had different GC and GC 3S contents indicating different codon usages. Among the zoonotic CoVs, SARS-CoV-2 had the lowest ENC and CAI values, suggesting that SARS-CoV-2 is poorly adapted to human codon usage. However, these values are similar to those of the less-pathogenic human-specific CoVs. In addition, we showed that the ENC values and the CAI values of SARS-CoV-2 were decreasing over first four months of the COVID-19 pandemic, which indicates that the codon usage bias of SARS-CoV-2 is increasing during the period, leading to further reduced pathogenicity.

The three zoonotic CoVs and the four human-specific CoVs appear to differ in codon usage. In this study, we compared the nucleotide compositions and the RSCU patterns among the eight groups of phylogenetically related CoVs. They all have AT-rich genomes, in agreement with the previous studies of CoVs (Gu et al., 2004; Jenkins and Holmes, 2003; Kandeel et al., 2020) . While, the genomes of human-specific CoVs have higher AT content than zoonotic CoVs, SARS-CoV-2 had the highest AT content among zoonotic CoVs. In correlation with the previous knowledge of CoV genome composition, CoVs frequently use A-and U-ended codons (Chen et al., 2017; Dilucca et al., 2020) . Human-specific CoVs have higher A/U frequency at the third codon position than zoonotic CoVs, indicating that the codon usage of zoonotic CoVs is different from the human-specific CoVs. In addition, SARS-CoV-2 has the highest A/U frequency at the third codon position among the three zoonotic CoVs. This indicates that the codon usage of SARS-CoV-2 is different from the other zoonotic CoVs. These were supported by the result of the PCA analysis of RSCU data, which placed these CoVs in different clusters in agreement with their phylogenetic relationships.

Among the three zoonotic CoVs, SARS-CoV-2 has the most codon preference. Although the ENC values of the zoonotic CoV and human-specific CoV genomes were all higher than 35, human-specific CoVs displayed higher codon usage biases than the zoonotic CoVs, and SARS-CoV-2 displayed a higher codon usage bias than the other zoonotic CoVs. In human-pathogenic viruses, high codon preference has generally been associated with reduced host range (Kumar et al., 2018) . Thus far, the four human-specific CoVs are known to infect only humans (Fehr and Perlman, 2015) . On the contrary, other than humans, SARS-CoV-2 can infect cats and ferrets . While the sequence characteristics of the receptor in the host is crucial in determining the susceptibility of animals to CoVs, the codon preference by SARS-CoV-2 could potentially make some animals inefficient hosts. Whether SARS-CoV-2 might have a narrower host range than the other two zoonotic CoVs remains to be determined.

SARS-CoV-2 appears to be poorly adapted to human codon usage compared with other zoonotic CoVs. Among the three zoonotic CoVs, SARS-CoV-2 had the lowest human codon CAI value, indicating that it has the lowest adaptation to the human host. Thus, it probably has the lowest likely efficiency of gene expression among the three zoonotic CoVs, as seen in comparative analyses of other viruses (Kumar et al., 2018; Zhang et al., 2019 ). This appears to contradict the finding in one J o u r n a l P r e -p r o o f Journal Pre-proof study (Dilucca et al., 2020) , but is in agreement with the finding in another (Kandeel et al., 2020) . On the other hand, the four human-specific CoVs have lower human codon CAI values than the zoonotic CoVs, and among the human-specific CoVs, the CAI values of HCoV-229E and HCoV-OC43 higher than that of HCoV-NL63 and HCoV-HKU1. These findings indicate that HCoV-229E and HCoV-OC43 have higher adaptation to the human host than the other two. Over the years, the prevalence of human-specific CoVs in adults and infants is in the following order: HCoV-OC43, HCoV-229E, HCoV-NL63, and HCoV-HKU1 (Bouvier et al., 2018; Zeng et al., 2018) . While these four CoVs cause only mild, self-limiting respiratory infections in humans, HCoV-OC43 infected patients have a significantly higher rate of acute respiratory tract infections than HCoV-HKU1 infected patients (Gerna et al., 2007) . Therefore, the outcome of the CAI analysis suggests that SARS-CoV-2 might have the lowest translation efficiency of viral proteins in humans among zoonotic CoVs, contributing to its lower pathogenicity than SARS-CoV and MERS-CoV (Chen, 2020; Wang et al., 2020) . In other RNA viruses, codon bias has been associated with reduced virulence (Diaz-San Segundo et al., 2016; Goñi et al., 2012; van Weringh et al., 2011) .

Mutational pressure probably plays a major role in shaping the codon usage in SARS-CoV-2. In ENC-GC 3S plots, most values clustered on or around the expected curves for all zoonotic and human-specific CoVs, which is an indication of the dominant role of mutational pressure in shaping codon usage of viruses (He et al., 2016) . The E gene of SARS-CoV-2, however, appears to be evolved differently from others. Compared with other zoonotic CoVs, the E gene of SARS-CoV-2 has the lowest ENC value. In all zoonotic CoVs, the E gene also has CAI values much lower CAI values to human codons than other genes. This suggests that the E gene of SARS-CoV-2 has a higher codon preference than other genes. Although mutational pressure plays a dominant role in shaping codon usage of all zoonotic CoVs, the E gene in SARS-CoV-2 appears to be under natural selection, as its ENC-GC 3S values were dispersed under the expected curve. The function of the E protein is not clear. It has been suggested to play an important role in viral assembly, budding, envelope formation, and pathogenesis (Schoeman and Fielding, 2019) .

Over the initial four months of the COVID-19 pandemic, SARS-CoV-2 appears to have increased codon preference. The ENC values of SARS-CoV-2 genomes collected over the 4-month period gradually decreased, indicating the codon usage preference has been increasing in SARS-CoV-2. This contradicts somewhat with the recent report indicating that nearly 80% of the recurrent mutations in SARS-CoV-2 produced non-synonymous changes at the protein level, possibly suggesting that adaptation of the virus to humans is occurring (van Dorp et al., 2020) . However, the reduced adaptation of the virus to human codons is collaborated in the present study by the decreasing CAI values over the 4-month period. The reduced ENC and CAI values are both suggesting that pathogenicity of SARS-CoV-2 in humans could be decreasing. In contrast, the E gene of the virus, which is the only gene under natural selection (indicated by ENC-GC 3S plot), seems to be adapting to human codons. Further studies of the E gene are needed to elucidate its role in the adaptative evolution of SARS-CoV-2.

In conclusion, results of the study suggest that the three zoonotic CoVs have different extent of adaptation to the human host. Among them, SARS-CoV-2 has the lowest ENC and CAI values, which are more similar to those from the human-specific Veit, M., Su, S., 2019. Genetic evolution and molecular selection of the HE gene of influenza C virus. Viruses 11. https://doi.org/10.3390/v11020167 Zhou, P., Yang, X. L., Wang, X. G., Hu, B., Zhang, L., Zhang, W., Si, H. R., Zhu, Y., Li, B., Huang, C. L., Chen, H. D., Chen, J., Luo, Y., Guo, H., Jiang, R. D., Liu, M. Q., Chen, Y ., Shen, X. R., Wang, X., Zheng, X. S., Zhao, K., Chen, Q. J., Deng, F., Liu, L. L., Yan, B., Zhan, F. X., Wang, Y. Y., Xiao, G. F., Shi, Z. L., 2020 . A pneumonia outbreak associated with a new coronavirus of probable bat origin. Nature 579, 270-273. https://doi.org/10.1038/s41586-020-2012-7 J o u r n a l P r e -p r o o f

Comprehensive analysis of genetic and evolutionary features of the hepatitis E virus

Causes and implications of codon usage bias in RNA viruses

Species-specific clinical characteristics of human coronavirus infection among otherwise healthy adolescents and adults

Modulation of poliovirus replicative fitness in HeLa cells by deoptimization of synonymous codon usage in the capsid region

Codon adaptation index as a measure of dominating codon bias

Dissimilation of synonymous codon usage bias in virus-host coevolution due to translational selection

Pathogenicity and transmissibility of 2019-nCoV-A quick overview and comparison with other emerging viruses

Analysis of the codon usage pattern in Middle East Respiratory Syndrome Coronavirus

An evaluation of measures of synonymous codon usage bias

Hepatitis A virus adaptation to cellular shutoff is driven by dynamic adjustments of codon usage and results in the selection of populations with altered capsids

Origin and evolution of pathogenic coronaviruses

SARS and MERS: recent insights into emerging coronaviruses

Synonymous deoptimization of foot-and-mouth disease virus causes attenuation in vivo while inducing a strong neutralizing antibody response

Codon usage and phenotypic divergences of SARS-CoV-2 genes

Coronaviruses: an overview of their replication and pathogenesis

Molecular evolution of human coronavirus genomes

Human respiratory coronavirus HKU1 versus other coronavirus infections in Italian hospitalised patients

Pandemic influenza A virus codon usage revisited: biases, adaptation and implications for vaccine strain development

Analysis of synonymous codon usage in SARS Coronavirus and other viruses in the Nidovirales

New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0

Analysis of codon usage patterns in Ginkgo biloba reveals codon usage tendency from A/U-ending to G/C-ending

The extent of codon usage bias in human RNA viruses and its evolutionary origin

From SARS and MERS CoVs to SARS-CoV-2: Moving toward more biased codon usage in viral structural and nonstructural genes

Analysis of Nipah virus codon usage and adaptation to hosts

Evolution of codon usage bias in henipaviruses is governed by natural selection and is host-specific

Identifying SARS-CoV-2-related coronaviruses in Malayan pangolins

Insights into the genetic and host adaptability of emerging porcine circovirus 3

Zika virus attenuation by codon pair deoptimization induces sterilizing immunity in mouse models

Estimates of the ongoing need for social distancing and control measures post-"lockdown" from trajectories of COVID-19 cases and mortality

COVID-19, SARS and MERS: are they closely related?

jModelTest: phylogenetic model averaging

Coronavirus envelope protein: current knowledge

An evolutionary perspective on synonymous codon usage in unicellular organisms

The codon Adaptation Index--a measure of directional synonymous codon usage bias, and its potential applications

Susceptibility of ferrets, cats, dogs, and other domesticated animals to SARS-coronavirus 2

Synonymous codon usage in Escherichia coli: selection for translational accuracy

Molecular epidemiological study on Infectious Pancreatic Necrosis Virus isolates from aquafarms in Scotland over three decades

Emergence of genomic diversity and recurrent mutations in SARS-CoV-2

HIV-1 modulates the tRNA pool to improve translation efficiency

Receptor recognition by the novel coronavirus from Wuhan: an analysis based on decade-long structural studies of SARS coronavirus

Clinical characteristics of 138 hospitalized patients with 2019 novel coronavirus-infected pneumonia in Wuhan, China

Codon usage bias and the evolution of influenza A viruses. Codon usage biases of influenza virus

Isolation of SARS-CoV-2-related coronavirus from Malayan pangolins

Epidemiology and clinical characteristics of human coronaviruses OC43, 229E, NL63, and HKU1: a study of hospitalized children with acute respiratory tract infection in Guangzhou, China