We obtained the 3-dimensional structures of full-length S and RBD-hACE2 of both SARS-CoV-1 and SARS-CoV-2 from the RCSB Protein Data Bank (PDB) (Protein Data Bank 2019). The structure of a trypsin-cleaved SARS-CoV-1’s spike glycoprotein (PDB ID: 6ACG) was used for stability analysis. The structure of the complex of SARS-CoV-1 S RBD and the human hACE2 receptor (PDB ID: 2AJF) was used for stability and interaction analysis. For SARS-CoV-2, we obtained the protein complex structure of RBD-hACE2 (PDB ID: 6LZG) for stability and interaction analysis. The SARS-CoV-2 S (PDB ID: 6VYB) was used only for stability analysis. For structural alignment of RBD of SARS-CoV-1 and SARS-CoV-2, we used the structures PDB ID:2AJF and PDB ID:6M17, respectively. For structural alignment of the S protein of SARS-CoV-1 and SARS-CoV-2, we used the structures PDB ID:6ACG and PDB ID:6VYB, respectively. PyMOL (Schrödinger, LLC 2015) was used for the visualization of the PDB structures and for structural alignments.

Foldx version 5 (Schymkowitz et al., 2005), coded by the Foldx Consortium, was used for mutational analysis. We used the command line interface of Foldx to mutate each residue to the other 19 residues. Foldx calculates free energy, ΔG, by using the contributions of hydrophobic and polar groups to the solvation energy, Van der Waals, hydrogen bonding, and electrostatic interactions. These energy parameters were experimentally derived (Schymkowitz et al., 2005). In this study, we used the default parameters for the computation of the wildtype and mutant free energies. The initial step used the ‘RepairPdb’ command to repair the wildtype protein structure. The ‘RepairPdb’ command finds the minimum energy conformation for the protein structure by flipping side chains of all residues, especially Asparagine, Glutamine, and Histidine to reduce steric clashes (Schymkowitz et al., 2005). This was followed by either the use of the ‘BuildModel’ command and the ‘AnalyseComplex’ command to calculate the folding energy change and the interaction or binding energy change, respectively. For each mutation, we used Foldx to calculate the folding energy change (ΔΔG) and binding energy change (ΔΔΔG) (FOLDX, 2021). The mathematical equation for the calculation of folding energy change (ΔΔG) is: Δ Δ G ( stability ) = Δ G ( folding ) MUT − Δ G ( folding ) WT

Theoretically, a negative ΔΔG means that the mutation leads to a more stable protein structure and a positive ΔΔG means that the mutation leads to a less stable protein structure. The five categories of the impact of the folding energy change (ΔΔG) are-highly stabilizing (ΔΔG < –2.0 kcal/mol), moderately stabilizing (–2.0 < ΔΔG < –0.5 kcal/mol), neutral (0.5 < ΔΔG < +0.5 kcal/mol), moderately destabilizing (+0.5 < ΔΔG < 2.0 kcal/mol), and highly destabilizing (ΔΔG >2.0 kcal/mol). The mathematical equation for the calculation of binding energy changes (ΔΔΔG) is: Δ Δ Δ G ( binding ) = Δ Δ G ( binding ) MUT − Δ Δ G ( binding ) WT

A negative ΔΔΔG means the mutation strengthens the binding energy and a positive ΔΔΔG means the mutation weakens the binding energy. The effect of the binding energy changes was also classified into five categories: large affinity decrease (ΔΔΔG >0.5 kcal/mol), moderate affinity decrease (0.1 < ΔΔΔG ≤0.5), neutral (–0.1 < ΔΔΔG ≤0.1 kcal/mol), moderate affinity increase (–0.5 < ΔΔΔG ≤ –0.5 kcal/mol), and large affinity in-crease (ΔΔΔG < –0.5 kcal/mol).

We used the Polymorphism Phenotyping v2 (PolyPhen2) (I. A. Adzhubei et al., 2010) and Screening for non-acceptable polymorphisms (SNAP) (Bromberg and Rost 2007) prediction tools to predict the pathogenicity of each missense mutations. We utilized the R programming language (https://www.r-project.org/) for data visualization for the purpose of drawing inferences. Specifically, we constructed boxplots to compare the prediction of pathogenicity between PolyPhen2 and SNAP.

The FASTA sequences of SARS-CoV-1 and SARS-CoV-2 S proteins were retrieved from the universal protein knowledgebase (UniProtKB) (Bateman 2019). We performed the pairwise sequence alignment of SARS-CoV-1 (Entry: P59594) and SARS-CoV-2 (Entry: P0DTC2) using the Clustal Omega computer program (https://www.ebi.ac.uk/Tools/msa/clustalo/) and Jalview2 (www.jalview.org). We performed structural alignments of the 3-D structures of SARS-CoV-1 and SARS-CoV-2 using PyMOL (http://www.pymol.org/). The “fetch” and the “align” commands on PyMOL aligned the single chains of the spike proteins of SARS-CoV-1 (PDB ID: 6ACG, chain A) and SARS-CoV-2 (PDB ID: 6VYB, chain A), and RBD of SARS-CoV-1 (PDB ID: 2AJF, chain E) and SARS-CoV-2 (PDB ID: 6M17, chain E).

We compared the outputs of six computational prediction tools on SARS-CoV-1 S stability. Each method utilizes different protein structural properties in predicting the effects of mutations on wild-type protein structures. The mutation Cut-off Scanning Matrix (mCSM) tool encodes atomic-distance patterns to predict the impact of mutations on protein structure (Pires et al., 2014a). The Site Directed Mutator (SDM) uses a statistical potential energy function to calculate a stability score. SDM applies a cut-off of 2 kcal/mol to classify stabilizing and destabilizing mutations (Pandurangan et al., 2017). DUET is a tool that combines, consolidates, and optimizes mCSM and SDM tools (Pires et al., 2014b). The overall accuracy of DUET is better than either mCSM or SDM (Pires et al., 2014b). The DynaMut tool utilizes Normal Mode Analysis (NMA) and graph-based signatures to predict the impact of missense mutations (Rodrigues et al., 2018). Finally, the I-mutant suite 3.0 tool uses a support vector machine (SVM) algorithm and accepts the protein sequence or structure as input. However, it predicts the impact of mutations more accurately with the protein structure inputted (Capriotti et al., 2005).

We also compared four computational prediction tools on SARS-CoV-1 S RBD affinity. The mutation Cut-off Scanning Matrix on Protein-Protein Interaction (mCSM-PPI2) tool uses inter-residue network complexes and graph-based signatures to predict the impact of missense mutations on protein affinity (Rodrigues et al., 2019). The Muta-Bind2 tool evaluates the changes caused by missense mutations on protein affinity. Muta-Bind2 is also a useful tool in protein design (Li et al., 2016). Lastly, PISA (Protein, Interfaces, Structures, and Assemblies) tool from the Protein Data Bank in Europe (PDBe) was used to analyze interface residues by comparing the contributions of their solvation energy to the interaction energy (Krissinel and Henrick 2005).

We modified the outputs from these tools to stay consistent with Foldx output. That is, significant positive ΔΔG and ΔΔΔG values destabilize and weaken binding affinity, respectively, and vice versa.

The Jalview tool (Waterhouse et al., 2009) shows the aligned residues, the quality of the alignment, the conservations scores, and the consensus between the RBD sequences of SARS-CoV-1 and SARS-CoV-2 (Figure 1A). The sequence alignment using Clustal Omega algorithm within the Jalview tool indicates a ∼76% sequence identity between RBD regions of SARS-CoV-1 and SARS-CoV-2. The bright yellow bars represent high conservation and quality between the aligned residues within each column. As shown in Figure 1B, the structural alignment revealed the evolutionary relationship between SARS-CoV-1 and SARS-CoV-2. SARS-CoV-1 S (PDB ID: 6ACG, chain A) aligned with SARS-CoV-2 S (PDB ID: 6VYB, chain A) with an RMSD of 2.272. Furthermore, we performed the structural alignment of the RBDs of SARS-CoV-1(PDB ID: 2AJF, chain E) and SARS-CoV-2 (PDB ID: 6M17, chain E), and it yielded an RMSD of 1.043. We selected a shorter fragment, RBM, from the RBD and performed the structural alignment of the RBM of both coronaviruses. This alignment yielded an RMSD of 0.878, which is more homologous. The shape and spatial orientation of the structural alignments overlapped which indicates close atomic coordinates between the two structures.

The SARS-CoV-1 S protein has 1,255 residues that were used to generate 23,845 non-redundant missense mutations. The effect of each mutation on the stability of the SARS-CoV-1 S protein was evaluated. Of the total mutations performed, 20,083 missense mutations generated energy changes, while the remaining 3,762 missense mutations gave no output due to missing residues on the protein structure. Figure 2A shows that 11,635 of the 20,083 (58%) missense mutations increased the free energy of the S protein by at least 0.5 kcal/mol, 2,964 of 20,083 (15%) missense mutations reduced the S protein’s free energy by at most –0.5 kcal/mol, and 5,484 of 20,083 (27%) had a neutral effect on the stability of the wildtype S protein. The standard error of energies calculation using the Foldx suite is ∼0.5 kcal/mol (Schymkowitz et al., 2005). Therefore, the folding energy changes within the range (–0.5 < ΔΔG <0.5) are insignificant or categorized as having neutral effect. In more specific categories, 4,767 mutations had a highly destabilizing effect (ΔΔG > 2.5 kcal/mol) on the spike protein, 6,868 mutations moderately destabilize the spike’s protein (0.5 < ΔΔG ≤ 2.5 kcal/mol), 5,484 mutations had a neutral effect (–0.5 < ΔΔG ≤ 0.5 kcal/mol), 2,816 mutations moderately stabilize the S protein (–2.5ΔΔG ≤ –0.5 kcal/mol), and 148 mutations have a highly destabilizing effect (ΔΔG < –2.5 kcal/mol) on the S protein.

The line chart and the overall heatmap show the distribution of the mutations along the entire length of the spike protein (Figure 3A). The red lines represent the positive mean values while the blue lines represent the negative mean values. The bubbles represent the folding energy changes when all the residues were mutated to Alanine. Based on the average of ΔΔG, mutations with the highest destabilizing effects were found in positions 839, 634, 418, 430, 536, and 1,113 (Figure 3). The two most destabilizing missense mutations, A430W and A430F, cause energy changes at 66.18 kcal/mol and 56.4 kcal/mol, respectively. This position is also within the RBM of the SARS-CoV-1 spike S protein. On the other hand, mutations with the highest stabilizing effects were found in positions: 1,059, 981, 500, 1,089, 150, and 247 (Figure 3). The two most stabilizing missense mutations, G981W and T1059F, reduced the free energy of the wildtype structure by –5.16 kcal/mol and –4.98 kcal/mol, respectively. The residue positions with the highest and lowest mean folding energy changes (∆∆G) were G839 and T1059, which are both within the S2 subunit of the S protein. Compared to the mean values of all mutations, the missense mutations to Alanine had more destabilizing effects than stabilizing effects. The two “white” gaps in the overall heatmap represent missing values in the residue positions 661–673 and 812–831, respectively.

We utilized the crystal structure of the RBD of SARS-CoV-1 (PDB ID:2AJF) for computational prediction of the effect of computed mutations on the stability of SARS-CoV-1 RBD. This analysis gave an alternative perspective on the predictive power of computational tools. We generated 3,841 mutations and calculated the folding energy change (∆∆G) caused by each missense mutation. Like the full spike analysis, 63% of the missense mutations destabilized the SARS-CoV-1 RBD structure. Meanwhile, 8% of the missense mutations stabilized the SARS-CoV-1 RBD structure. Figure 2B shows the pie chart of the effects of the missense mutations on stability of SARS-CoV-1 S RBD. The distribution of the effect of all missense mutations in both the RBD and the full-length S stability analysis correlates (r = 0.6824). Observation of individual mutations revealed analytic similarity in both stability analysis. The heatmaps in Figures 3B, 4A also show similarity in mis-sense mutations with large effects on the stability of the SARS-CoV-1 S stability. For example, the missense mutation A430W had a large effect on both the stability of the RBD (Figure 4B) and the entire spike protein of the SARS-CoV-1.

We compared the mutational effects of seven residues on SARS-CoV-1 S protein to corresponding residues on SARS-CoV-2 S protein. Four residues (T1059, G981, S500, and R1089) had the highest mean stabilizing effects, while three residues (G839, G634, and A430) had the highest mean destabilizing effects. Pairwise alignment of SARS-CoV-1 S and SARS-CoV-2 S showed that the residues T1059, G981, S500, R1089, G839, G634, and A430 on SARS-CoV-1 corresponds with residues T1077, G999, S515, R1107, G857, G648, and S443 on SARS-CoV-2, respectively. As shown in Table 1, the missense mutations on SARS-CoV-1 S protein had similar effects on its stability when compared to corresponding residues on SARS-CoV-2 S protein. We observed that the missense mutation A430S destabilizes SARS-CoV-1 S protein by introducing the folding energy change at 1.8 kcal/mol. Consistently, S443A in corresponding position of SARS-CoV-2 increase the S stability (∆∆G = –0.768 kcal/mol). Interestingly, SARS-CoV-1 D600G can destabilize S (∆∆G = 0.21 kcal/mol). D614G, the dominant variant of SARS-CoV-2, corresponds to D600G but was predicted to stabilize S (∆∆G = –0.784 kcal/mol). Another In position 247, the substitution of Threonine with Alanine increased the stability of the SARS-CoV-1. However, the corresponding residue on SARS-CoV-2 A260 was missing. Furthermore, we compared the predicted stability effects of target SARS-CoV-2 S RBD mutations with the current deep mutational scanning approach (Starr et al., 2020). The computational predictions agree with the experimental results (Table 1).

The protein structure (PDB ID:2AJF) used for interaction analysis covered the S RBD that interact with the hACE2 receptor in humans. This RBD chain contains 180 residues, which spanned from residual position 323 to 502. A total of 3,420 mutations were computed, with 114 missing values. The remaining 3,306 mutations were classified into one of five categories according to their binding energy changes (ΔΔΔG). As shown in Figure 2C, 257 mutations had a large-decrease effect, 183 mutations had a moderate-decrease effect, 2,738 mutations had a neutral effect, 97 mutations had a moderate-increase effect, and 31 mutations had a large-increase effect. Figure 5E shows the mutations with the largest effects, based on ΔΔΔG mean values from all possible mutations in that position.

The seven topmost mutations with decreasing effects are found in positions which have mean values greater than 1.5 kcal/mol. While the six topmost mutations with in-creasing effects are found in positions which have mean values less than –0.3 kcal/mol. The values in red boxes represent the minimum values in each position, while the values in blue boxes represent maximum values. The overall distribution of the missense mutations as they change ACE2 binding affinity of the SARS-CoV-1 S protein can be seen in Figure 6. There are 7 distinct regions or domains in the heatmap: [337–338], [388–396], [402–409], [424–436], [439–447], [460–467], and [472–494]. Four of the seven distinct regions are concentrated within the RBM, and they span longer stretches of residues. Most of the red and blue spikes on the line chart are also within the RBM. These distinct regions can provide insight into which residues play a key role in the interaction between SARS-CoV-1 S protein and the hACE2 receptors in humans.

As shown in Table 1, Asparagine is found on position 501 in SARS-CoV-2 which corresponds to Threonine on position 487 in SARS-CoV-1. The computed missense mutation T487N increased the ∆∆∆G by 2.28 kcal/mol and decreased the ∆∆G by –0.14. In contrast, N501T, on the SARS-CoV-2 decreased ∆∆∆G by –1.346 kcal/mol. Tyrosine and Tryptophan have nonpolar aromatic chains therefore the missense mutations T487Y and T487W, had similar effects (decrease binding affinity) on corresponding residue on SARS-CoV-2, N501Y and N501W, respectively. Most of the missense mutations on positions S432, P462, and N479, on the average, increase the binding affinity of SARS-CoV-1. However, the corresponding residues on SARS-CoV-2 are Valine, Alanine, and Glutamine, respectively. The effect of missense mutations on these residues were evaluated by examining missense mutations on corresponding positions on SARS-CoV-1. However, mis-sense mutation N479Q decreased the ∆∆∆G by –0.607 kcal/mol and increased ∆∆G by 1.91 kcal/mol. A possible explanation is that the missense mutations created Van der Waal clashes with neighboring atoms because of shorter distance; The missense mutation, T487Y, caused a differential Van der Waals of -0.449 kcal/mol, which led to a greater re-pulsion from nearby atoms. Recent SARS-CoV-2 variant L452R (∆∆∆G = –0.395 kcal/mol; ∆∆G = 0.021 kcal/mol) corresponds to SARS-CoV-1 K439R (∆∆∆G = 0.247 kcal/mol; ∆∆G = 0.41 kcal/mol). The change in residue from Lysine in SARS-CoV-1 S protein to Leucine in SARS-CoV-2 S protein may be responsible for the increase in binding affinity caused by L452R.

Post-Translational modifications (PTMs) of the SARS-CoV-1 S protein are responsible for the folding, maturation, and function of the S protein. An important PTM utilized by SARS-CoV-1 is the O- and N-linked Glycosylation, which plays a key role in the shielding of viruses from the host’s immune system (Watanabe et al., 2019). As a result, viruses evolve to become glycosylated as much as possible (Sugrue 2007). One O-linked glycosylation site, S336, located within the RBD was predicted using online webserver (Steentoft et al., 2013). On the average, the computed missense mutations in this position increase the stability of SARS-CoV-1. For example, the missense mutation, S336E, reduces the ∆∆G by -1.93 kcal/mol, thereby stabilizing the SARS-CoV-1 spike protein. Conversely, S336V in-creases the ∆∆G by 1.28 kcal/mol, thereby destabilizing the SARS-CoV-1 structure. We used five online prediction tools in our N-linked glycosylation sites prediction- NetNGly 1.0 (Kumar et al., 2020), N-GlyDE (Pitti et al., 2019), SPRINT-Gly (Taherzadeh et al., 2019), Glycopp v1.0 (Chauhan et al., 2012), and Glycopred (Hamby and Hirst 2008). We narrowed our results to sites predicted by at least two or three prediction tools. Of these, we selected 20 N-linked glycosylation sites located within the spike protein of the SARS-CoV-1. Three putative N-linked glycosylation sites, N318, N330, N357, were located within the RBD, and they were predicted to have no effect on the binding affinity of SARS-CoV-1 spike protein to hACE2. On the average, mutations at position N318 were predicted to increase the stability of SARS-CoV-1. However, on the average, the mutations on N330 and N357 were predicted to have neutral effects.

Figure 7 shows the target mutations in twenty N-linked glycosylation sites, one O-linked glycosylation site, and one palmitoylation site. Palmitoylation modifies the spike protein through cysteine-rich residues (Petit et al., 2007). This is known to mediate the fusion of the spike protein to hACE2. Predicted palmitoylation sites (Ren et al., 2008) were found in positions C19, C1217, C1218, C1222, C1232, C1235, and C1236. The average computed mutations at the C19 position predicted neutral effects, which showed no significant effect in the stability of SARS-CoV-1. However, the remaining six palmitoylated were located outside the residues covered by the structural protein.

We used sequence-based mutation pathogenicity tools to predict the damaging effect of our computed mutations on the SARS-CoV-1 S function. We analyzed 25,101 missense mutations using the full-length SARS-CoV-1 S (1–1,255). The PolyPhen2 scores gave probabilistic values on the tolerance and deleterious effect of missense mutations. A score less than 0.446 is considered benign, a score greater than 0.446 but less than 0.908 is considered possibly damaging, and a score greater than 0.908 is considered probably damaging. In Figure 8A, the missense mutations with neutral effect were predicted to be mostly tolerated with some classified as benign, while the mean value, as shown by the red line, is considered possibly damaging. Whereas the mean values of the moderately increasing and decreasing mutations were predicted to be possibly damaging. However, the mean values of large increasing and decreasing mutations were predicted to be probably damaging. The analysis of variance (ANOVA) showed that the means of all five categories were significantly different, with p-value < 2e-16 (Figure 8A).

SNAP is a neural network machine learning algorithm that accepts protein structure as input for functional predictions (Bromberg and Rost 2007). SNAP scores less than zero have neutral effect, while SNAP scores greater than zero have a pathogenic effect. 13,179 (∼53%) of the total 25,101missense mutations have pathogenic effects. We performed further analysis to correlate the effects of folding energy change ∆∆G on stability with their SNAP scores. As shown in Figure 8B, the SNAP scores of missense mutations with effects greater than 2.5 kcal/mol or less than –2.5 kcal/mol were higher than missense mutations with moderate effects (0.5 < ΔΔG ≤ 2.5 or –2.5 =< ΔΔG < –0.5 kcal/mol). Furthermore, the missense mutations with neutral effects on SARS-CoV-1 stability also gave a neutral SNAP prediction. The statistical analysis showed that the correlation of SNAP scores in the five groups were significant with a p-value of <2e-16.

This study analyzed mutations, found on the SARS-CoV-1 protein, which have been verified through experiments to have cellular or molecular effects on its functions. Analysis of SARS-CoV-1 isolates from the 2002–2003 and 2003–2004 outbreaks revealed the main driving mutations (Li W. et al., 2005). Six residue positions of S protein, 344, 360, 472, 479, 480, and 487, were highlighted to have varying residues between the two isolates (Li W. et al., 2005). The 2002–2003 isolate had residues K344, F360, L472, N479, D480, and T487. Meanwhile, the 2003–2004 isolate had the following corresponding residues R344, S360, P472, K479, G480, and S487. The effects of the residue changes from the 2002–2003 isolate to the 2003–2004 isolate are shown in Table 2. The residue change, N479K, would increase the binding affinity (∆∆∆G = –1.008 kcal/mol) and reduced S stability (∆∆G = 1.52 kcal/mol). Two mutations, L472P and T487S, decreased the binding affinity of S RBD-ACE2 by introducing ∆∆∆G at 1.875 kcal/mol and 0.752 kcal/mol, respectively. Furthermore, the mutation L472P destabilizes SARS-CoV-1 (∆∆G = 1 kcal/mol). Therefore, understanding how the changes in these residues affect the stability and infectivity of the SARS-CoV-1 in the two isolates would help in the target of specific residues on SARS-CoV-2. Recent research to study the mutation of residues from the 2002–2003 isolates to residues in 2003–2004 isolate is consistent with this computational study (Li W. et al., 2005).