Sequence analysis of SARS-CoV-2 genome reveals features important for vaccine design

As the SARS-CoV-2 pandemic is rapidly progressing, the need for the development of an effective vaccine is critical. A promising approach for vaccine development is to generate, through codon pair deoptimization, an attenuated virus. This approach carries the advantage that it only requires limited knowledge specific to the virus in question, other than its genome sequence. Therefore, it is well suited for emerging viruses, for which we may not have extensive data. We performed comprehensive in silico analyses of several features of SARS-CoV-2 genomic sequence (e.g., codon usage, codon pair usage, dinucleotide/junction dinucleotide usage, RNA structure around the frameshift region) in comparison with other members of the coronaviridae family of viruses, the overall human genome, and the transcriptome of specific human tissues such as lung, which are primarily targeted by the virus. Our analysis identified the spike (S) and nucleocapsid (N) proteins as promising targets for deoptimization and suggests a roadmap for SARS-CoV-2 vaccine development, which can be generalizable to other viruses.


Scientific Reports
| (2020) 10:15643 | https://doi.org/10.1038/s41598-020-72533-2 www.nature.com/scientificreports/ the codon usage of their hosts 12,13 , a phenomenon that is not well understood. In this regard, a thorough characterization of codon, codon pair and dinucleotide usage of SARS-CoV-2 can provide useful information regarding expression potential of the viral genes and the fitness of the virus in its human or other hosts. Furthermore, it has been shown that viral attenuation can be achieved through extensive changes in codon pair usage of viral genes 2 . Since the mechanism of viral attenuation through codon pair deoptimization is not entirely clear, this in-depth analysis is necessary to guide the development of new vaccines. Coronaviruses (CoVs) are enveloped, positive-stranded RNA viruses with a large genome of about 30 kb encoding multiple proteins 14 . Translation of a positive-stranded RNA from the initial infectious virus particles generates (among other proteins) a virally encoded RNA-dependent RNA polymerase (replicase). This replicase is necessary for viral replication and subsequent generation of viral subgenomic RNAs (sgRNAs), from which the synthesis of structural and accessory proteins occurs 14 . ORF1ab, which encodes the replicase polyprotein (among other proteins) occupies about two thirds of the 5′ prime end of this genome 14,15 . A − 1 programmed ribosomal frameshift (PRF) occurs half-way through ORF1ab, allowing the translation of ORF1b 14 . The efficiency of the frameshift thus modulates the relative ratios of proteins encoded by ORF1b and the upstream ORF1a and is critical for coronavirus propagation. Frameshift efficiency (ranging from 15 to 60%) in − 1 PRFs is commonly regulated by pseudoknotted mRNA structures following the frameshift, and the conservation of a three-stem pseudoknot in coronaviruses has been previously characterized 16 . Following ORF1ab, are the spike (S), ORF3a, envelope (E), membrane (M), ORF6, ORF7a, ORF7b, ORF8, nucleocapsid (N) and ORF10 genes. The S protein promotes attachment and fusion to the host cell, during infection 17 . In the case of SARS-CoV-2, S binds to the human angiotensin-converting enzyme 2 (ACE2) 15,18,19 . The E protein is an ion channel and regulates virion assembly 20 . The M protein also participates in virus assembly and in the biosynthesis of new virus particles 21 , while the N protein forms the ribonucleoprotein complex with the virus RNA 22 and has several functions, such as enhancing transcription of the viral genome and interacting with the viral membrane protein during virion assembly 23 . Many of the other ORFs have unknown functions or are not well characterized 24 , as their presence is not consistent across all coronaviruses.
We have conducted a thorough analysis of the codon, codon pair and dinucleotide usage of the SARS-CoV-2 and have assessed how it relates to other coronaviruses, its hosts, and to the tissues that SARS-CoV-2 has been reported to infect [25][26][27] . We have taken advantage of our recently published databases, which include genomic codon usage statistics for all species with available sequence data, and transcriptomic codon usage statistics from several human tissues 4,5,28 . We further analyzed each viral gene in terms of its codon characteristics and used an array of codon usage metrics that informed us of the potential of each gene sequence to contribute to the deoptimization of the virus. In the case of ORF1ab, we further examined the structure of the mRNA in the region following the frameshift, finding the SARS-CoV-2 mRNA to exhibit a similar pseudoknotted structure to known coronaviruses. We identified two viral genes, S and N, that represent valuable targets for deoptimization to generate an attenuated virus. In the future we plan to continue this in silico study with an experimental investigation of deoptimized virus and test whether the deoptimized S and/or N proteins are immunogenic, produce epitopes that are neutralizing, and result in antibodies that are escape-resistant. We believe, that our combined analysis can be used as a pipeline to guide codon pair deoptimization for viral attenuation and vaccine development or a posteriori to evaluate the effectiveness of the attenuation of a viral sequence.

Results
SARS-CoV-2 proximity to coronaviruses, host genomes and tissue transcriptomes. Since the end of last year when it first emerged, SARS-CoV-2 has been mutating and spreading around the world. Over 5,000 complete or near-complete SARS-CoV-2 genomes are currently accessible in GenBank, with various mutations. To determine which SARS-CoV-2 sequence was most appropriate to use, we retrieved all the published sequences of the virus available in NCBI's SARS-CoV-2 data hub (5,064 complete SARS-CoV-2 genomes); after excluding incomplete and low-quality sequences and CDSs with insertions or deletions, we calculated the percent difference in codon usage between these and the reference sequence. The average percent difference in codon usage was ~ 0.08%, or ~ 8 codons/10,000, clearly showing that variation in sequences is not significantly affecting overall codon usage. This degree of mutation between strains is corroborated by a recently published study 29 , and is encouraging as it suggests that escape mutants are unlikely to develop, even for viral genes that have the highest selection pressure such as the S protein. Furthermore, we examined genetic diversity data from Nextstrain 30 , accounting for 4,675 SARS-CoV-2 genomes. From these sequences, there are 293 codon positions (~ 23% of the S gene) with reported Shannon entropies, i.e. with a documented mutation at that position. Among all other genes, there are 2,328 such positions (~ 27% of all non-S genes), indicating a higher percentage of mutated codons outside of the S gene.
A discrepancy between virus and host codon and codon pair usage bias has been observed across a range of viruses 12,13,31,32 , therefore we examined whether this was true for SARS-CoV-2 and its current host and to other Coronaviruses. Codon pair data inherently contain the codon usage data and therefore are better suited than codon usage data for this type of comparison. As expected, SARS-CoV-2 codon pair usage closely resembles the codon usage of the coronoviridae family, while it is quite distinct from the codon pair usage of the human genome (Table 1). Bat (Chiroptera) and pangolin (Pholidota) from which the virus may have been transmitted to humans, as well as dog (Canis lupus familiaris) to which the virus is feared may be transmitted next, were included in the analysis. We find that these species have a similar codon usage when compared with human; therefore, viral tropism cannot be inferred based on codon usage data alone (Table 1). Since SARS-CoV-2 infects bronchial epithelial cells and type II pneumocytes and our recent findings show that transcriptomic tissue-specific codon pair usage can vary greatly from genomic codon pair usage 28 , we also examined the transcriptomic codon pair usage of the lung and how it compares with the SARS-CoV-2 codon pair usage. Rather surprisingly, the codon Scientific Reports | (2020) 10:15643 | https://doi.org/10.1038/s41598-020-72533-2 www.nature.com/scientificreports/ pair usage in the lung was more distinct from SARS-CoV-2 codon pair usage than the Homo sapiens genomic codon pair usage. The transcriptomic codon pair usage of kidney and small intestine, tissues that are also susceptible to the infection, are similarly distant from SARS-CoV-2 (Table 1). Recently, it was argued that some degree of dissimilation in codon usage between the virus and the host may be beneficial to the virus, as it does not severely impede host gene translation 11 .

Codon, codon pair and dinucleotide usage of SARS-CoV-2.
To inspect the sequence features of SARS-CoV-2 in more detail, we plotted its codon usage per amino acid and compared it with the human genome and lung transcriptome (Fig. 1). SARS-CoV-2 clearly exhibits a preference in codons ending in T and A (71.7%), which is not observed in the human genome (44.9% ending in T or A) and lung transcriptome (37.6% ending in T or A). Similarly, the kidney and small intestine transcriptome show a preference for codons ending in C and G (62.5% in the kidney and 61.8% in the small intestine, Supplemental Figure 1). The codon pair usage of SARS-CoV-2 was also examined in juxtaposition with the human codon pair usage ( Fig. 2A,B). The differences in codon pair usage of the two genomes are highlighted in Fig. 2C.
Since the mechanism of viral attenuation through codon pair deoptimization is not entirely clear, and it has been argued that it is an indirect result of increased CpG content, we further investigated the dinucleotide and junction dinucleotide profile of the SARS-CoV-2 as it compares with Homo sapiens genome and lung transcriptome (Fig. 3). Clearly, CpG dinucleotides are avoided in the SARS-CoV-2 genome, and to a lesser extent CC and GG dinucleotides are too. This provides an opportunity to increase immunogenicity of a potential attenuated virus vaccine by increasing its CpG content. RNA folding. The genome sequence determines not only the amino acid sequence, but also the structure of the mRNA. The mRNA structure following the frameshift site is expected to be especially biologically relevant, as pseudoknots following programmed ribosomal frameshifts have been found to regulate the efficiency of the frameshift 33,34 . We therefore sought to study the similarity of the SARS-CoV-2 mRNA structure compared with the structures of different coronavirus mRNAs in the region following the ORF1ab frameshift. Table 1. Euclidean distances between codon pair usage frequencies. Euclidean distance (scaled/1,000) between codon pair usage frequencies of SARS-CoV-2, Coronaviridae (All CoV), Homo sapiens (genomic), lung, kidney (cortex), small intestine (terminal ileum), Pholidota (pangolins), Chiroptera (bats) and Canis lupus familiaris. www.nature.com/scientificreports/ RNA structures were predicted using two distinct secondary structure prediction algorithms, LandscapeFold 36 and NuPack 35,37 . Of the top 10 coronaviruses whose predicted minimum free energy (MFE) structures best aligned to that of SARS-CoV-2, seven matched among the two algorithms, showing a high degree of agreement www.nature.com/scientificreports/ among the two sets of structure predictions. Those seven consensus best-aligned structures are shown, alongside the novel coronavirus post-frameshift structure, in Fig. 4A-H. The similarity of two of these structures to SARS-CoV-2 can be explained by a high degree of sequence similarity to the SARS-CoV-2 mRNA (a SARS-related coronavirus and a bat coronavirus, shown in Fig. 4B, C). However, the other five-all belonging to avian coronaviruses, which are part of the group of the so-called gammacoronaviruses, causing highly contagious diseases of chickens, turkey and other birds-were not in the top 10 sequences most closely aligned to the SARS-CoV-2 mRNA on the basis of sequence. It should be noted that, of the 5,064 SARS-CoV-2 sequences analyzed, 3,978 had a complete ORF1ab with the exact 'UUU AAA C' frameshift sequence in the annotated position. Of these, 3,951 share the same sequence in the 100 nts downstream of the signal, indicating a high degree of conservation in this region. Finally, we used LandscapeFold 36 to study the RNA folding beyond the MFE structures. We find that even those coronaviruses whose MFE structure does not contain a pseudoknot will fold into a pseudoknot in a relatively high fraction of cases, and that most coronaviruses have a relatively high probability of the initial stem following the frameshift folding into part of a 3-stem pseudoknot like the one exhibited by the SARS-CoV-2 MFE structure (Fig. 4I). Viral gene codon usage properties. We next sought to examine each viral gene separately in terms of their codon and codon pair usage. Relative synonymous codon usage (RSCU) and codon pair score (CPS) are commonly used metrics to describe the codon and codon pair usage bias, respectively. RSCU expresses the observed over expected synonymous codon usage ratio, while CPS is the natural log of the observed over expected synonymous codon pair ratio using observed individual codon usage 2,38 . In our analyses, RSCU and CPS are derived from human genomic codon and codon pair usage frequencies. For ease of comparison, we used ln(RSCU) to measure the codon usage bias. The average CPS across a gene is referred to as codon pair bias (CPB) of the gene 2 . The average ln(RSCU) and CPB of each viral gene was calculated and compared with host genes average ln(RSCU) and CPB (Fig. 5). The average RSCU, ln(RSCU) and CPB of each viral gene appear in Table 2. ORF10 was strikingly the least similar gene to the human genome in terms of both its codon and codon pair usage, followed by the E gene. These genes provide little opportunity for deoptimization, since their sequence is already far from optimal. On the other hand, genes S and N are more similar to human in terms

Discussion
We performed a comprehensive characterization of the codon, codon pair and dinucleotide usage of the SARS-CoV-2 genome with the intention to identify the best targets for codon pair deoptimization in order to design an attenuated virus for vaccine development. Genes N and S were singled out as the best potential targets for deoptimization, due to the relatively high CPB among the viral genes. Furthermore, they are structural proteins with known functions, which will facilitate subsequent studies. Currently, most published attempts of viral attenuation through codon pair deoptimization do not discuss the strategy for selecting which genes to deoptimize. Although codon pair deoptimization has been proven successful for viral attenuation 2,39,40 , the mechanism is not clear. In addition, it is reasonable to assume that for every successful published deoptimization attempt, there may be several unsuccessful and therefore unpublished ones. Similarly, there have been successful attempts to generate attenuated viruses through codon deoptimization [41][42][43][44] . Understanding the mechanism that leads to viral attenuation requires a thorough characterization of the viral sequence and of the consequences of sequence changes. There are several factors that may contribute to the efficacy of deoptimization strategies. In changing the codon pair usage, the dinucleotide frequency and the GC content are altered; mRNA secondary structure and translational kinetics are also perturbed. Further, the CpG content is changing, leading potentially to altered immunogenicity. It is likely that codon pair (or codon) deoptimization leads to reduced expression, either due to changes in transcription, mRNA stability, or translation efficiency 45,46 . Alternatively, it is possible that deoptimization may lead to perturbed cotranslational folding 47 , resulting in altered protein conformation. In the case of the S protein, this may lead to decreased binding affinity for the ACE2 protein, thus affecting viral fitness.  www.nature.com/scientificreports/ A number of parameters were considered in determining which proteins could be targets for codon pair deoptimization. It has been shown that deoptimizing one third or less of the virus is sufficient to attenuate the virus, and more extensive deoptimization may lead to a completely inactive virus 48 . ORF1ab takes up about two thirds of the virus; therefore, its size may make it an unsuitable target. Furthermore, altering its RNA sequence is likely to disturb the pseudoknot that is responsible for frameshifting. Since we have identified the sequence that is responsible for the frameshift, a partial codon pair deoptimization is possible. However, ORF1ab is essential for genome replication, which also does not support its capacity as a codon deoptimization target.
ORF10 has strikingly low CPB and RSCU; given that it is at the very end of the viral genome, there may be a structural reason for its nucleotide sequence. The E gene also has a very low CPB; interestingly, although ORF10 has both positive and negative CPS across its sequence, E has mostly negative CPSs, which make it unsuited for codon pair deoptimization. ORF7a is unusual, as it has the highest RSCU of the viral genes, but a rather low CPB (it uses preferred codons in unusual combinations). Although ORF7a is not the most compelling target for codon pair deoptimization, if a codon deoptimization strategy is attempted, this gene should be considered. It should, however, be noted that since it overlaps for a few nucleotides with ORF7b, any sequence changes should be considered in coordination for both genes.
Our sequence analysis pointed to the S and N genes as potential targets for codon pair deoptimization. The S protein, which binds to the cell-surface receptor and induces virus-cell membrane fusion, has the highest CPB score of all viral proteins, leading to significant flexibility for codon pair deoptimization. It should be noted that although protein S is a surface protein and is expected to be under pressure to avoid the immune system, an examination of the variants that have emerged over the past few months indicates that the S gene does not appear to mutate at a faster rate than other SARS-CoV-2 genes. The N protein, which forms the ribonucleoprotein complex with the virus RNA, is the most conserved and stable protein among the coronavirus structural proteins. It uses mostly codon pairs with intermediate frequency; thus, it could be substantially codon pair deoptimized.
The next step, which is beyond the scope of this in silico study, is to construct the deoptimized virus and test its infectivity, ability to replicate, and whether it retains any pathogenicity. The study should investigate whether deoptimized S and/or N proteins are immunogenic, produce epitopes that are neutralizing, and result in antibodies that are escape-resistant. While testing the fitness of the virus, our strategy of selecting targets would be tested and validated, which could lead to better understanding of the factors that make codon pair deoptimization successful in generating attenuated viruses. While these planned steps are essential in the experimental validation of our strategy, we aimed to provide the scientific community with the framework that would allow anyone to independently proceed towards engineering a codon deoptimized SARS-CoV-2 for vaccine development.
The risk of a new emerging virus is always present, and this has been poignantly highlighted by the current SARS-CoV-2. The current work could be used for the quick generation of a SARS-CoV-2 vaccine but also as a pipeline to facilitate vaccine development when the next virus is presented.

Materials and methods
Sequence accession and codon comparison. The complete reference sequences for severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2 , accession NC_045512.2) were downloaded from NCBI RefSeq 49 on June 8, 2020. CDS sequences from 5,064 complete SARS-CoV-2 isolates were downloaded using NCBI SARS-CoV-2 data hub on June 8, 2020 50 . Sequences of poor quality or with CDS lengths that did not match those of the reference sequence due to deletion or insertion were removed. To calculate percent difference in codon usage, each CDS was compared at the codon level to that of the reference sequence. Codons containing nucleotides where a base call could not be made ("N") were removed from the calculation. All scripts for this calculation were written in Python 3.7.4.
Comparison of codon and codon pair usage in host species. Codon, codon pair and dinucleotide usage data for Homo sapiens, Canis lupus familiaris, Chiroptera (bats) and Pholidota (pangolins) were downloaded from the CoCoPUTs database 5 on March 13, 2020. Likewise, human lung, kidney (cortex) and small intestine (terminal ileum) tissue-specific codon, codon pair and dinucleotide usage data were accessed from the TissueCoCoPUTs database 28 on March 13, 2020. Codon, codon pair and dinucleotide usage data for SARS-CoV-2 was calculated from the reference sequence (accession NC_045512.2) using scripts written in Python 3.7.4. Euclidean distances between codon pair usage frequencies were calculated using the dist function from the stats package in R 3.6.1. RNA folding. To ensure our results are robust to the prediction algorithm chosen as well as to the size of the window examined, we used two secondary structure prediction algorithms on two window sizes. We used NuPack to predict the minimum free energy (MFE) secondary structure on the 100 nucleotides (nts) following the frameshift 35 , and our own recently published free energy landscape enumeration algorithm, LandscapeFold, to examine the full structure landscape of the 75 nts following the frameshift 36 . For the latter, we employed the heuristic that the minimum stem length was set to 4. Aside from this heuristic, the two algorithms differ primarily in the loop entropy calculation, which especially affects the probability of pseudoknot formation.
Structure alignment was measured using a method similar to our previously-studied "per-base topology" score 36 . Taking the dot-bracket representation of each secondary structure, we summed the number of positions containing identical elements. Employing an alignment model allowing for gaps, with a gap penalty and a misalignment penalty of − 1, did not change our results 51 .