SARS-CoV-2 Spike mutations modify the interaction between virus Spike and human ACE2 receptors.

The high mutability of the SARS-CoV-2 virus is a growing concern among scientific communities and health professionals since it brings the effectiveness of repurposed drugs and vaccines for COVID-19 into question. Although the mutational investigation of the Spike protein of the SARS-CoV-2 virus has been confirmed by many different researchers, there is no thorough investigation carried out at the interacting region to reveal the mutational status and its associated severity. All the energetically favorable mutations and their detailed analytical features that could impact the infection severity of the SARS-CoV-2 virus need to be identified. Therefore, we have thoroughly investigated the most important site of the SARS-CoV-2 virus, which is the interface region (Residue 417-505) of the virus Spike that interacts with the human ACE2 receptor. Further, we have utilized molecular dynamic simulation to observe the relative stability of the Spike protein with partner ACE2, as a consequence of these mutations. In our study, we have identified 52 energetically favorable Spike mutations at the interface while binding to ACE2, of which only 36 significantly enhance the stabilization of the Spike-ACE2 complex. The stability order and molecular interactions of these mutations were also identified. The highest stabilizing mutation V503D confirmed in our study is also known for neutralization resistance.

Introduction

Viruses including SARS-CoV-2 change themselves in due course of time to sustain their existence; however, while most of the changes have very little or no effect at all, a few can severely impact the virus' epidemics including transmission pathogenicity, immune resistance, vaccine performance, and diagnosis. Therefore, it is of utmost importance to understand the virus’ organization and the mutational impact on it [1].

The genomic organization of SARS-CoV-2 consists of single-stranded positive-sense RNA of ∼30,000 bases that encode several proteins including structural, non-structural, and accessory [2,3]. Structural proteins, Spike (S), function for virus attachment and entry into the host. Envelope small membrane protein (E) is associated with virus budding; Membrane protein's (M) functional role is in virion morphogenesis, while viral genome packaging is done by nucleoprotein (N). Non-structural protein Nsp1-16 is obtained from precursor protein pp1a and pp1b, and their functional role is in viral replication, transcription, and production of enveloping proteins. Accessory proteins of SARS-CoV-2 (ORF3a, b, c, d, 6, 7a, b, 8, 9b, c, and 10) are responsible for the differential pathogenesis in a wide host range [4].

So far, several mutations in the proteins of SARS-CoV-2 have been identified that disrupt the normal functionality of the virus. For instance, in a replicating virus, the M protein is thought to be involved in enhanced glucose uptake and is therefore associated with virus replication, proliferation, and immune evasion [5]; any mutation in the M protein might change the rate of replication, proliferation, and immune evasion processes. The M protein also contains protein kinase C (PKC) sites, and mutation in the M protein could change these PKC sites, therefore impacting the dynamic ruffling and endocytosis [6]. Mutation in the E protein is associated with either enhancement or diminution of the viral pathogenesis persuaded via separate mechanisms. The former occurs when the mutation in the E protein favoring the replication and propagation and tight binding of the E protein to tight junction PLAS1 leads to inflammation, thus increasing pathogenesis [7]. The latter happens when mutation downgrading the E channel activity reduces the production of SARS-CoV-2 particles, therefore causing host cell death [8].

Spike mutation N439K prompted the emergence of the first lineage of SARS-CoV-2 B.1.141 and subsequently prompted the lineage B.1.258 [9,10]. Another Spike mutation Y453F enhances the Spike-ACE2 interaction found to be associated with lineage B.1.1.298, which also harbors the deletion of amino acid (AA) residues 69 to 70 at the N-terminal domain, emerges several times across the globe, and is highly infectious [11]. Serological analysis has determined that ∼90% of plasma or serum neutralizing antibodies target the RBD of the Spike, and therefore the reason for this immunodominance of RBD of the Spike could be a relative lack of glycan shielding [12,13].

In addition to immunodominance, escape mutation resisting neutralization is mainly located in the RBD domain of Spike identified by exposing the circulating variants and amino acid (AA) substitution variants to monoclonal antibodies and convalescent plasma containing polyclonal antibodies [11]. In consideration of the aforementioned facts, Spike mutations are very important and are hence receiving considerable attention, and several findings have been put forth, such as the idea that mutation in the Spike protein can cause several consequences such as stabilization/destabilization in virus attachment to host ACE2, increased resistance to neutralizing antibody and convalescent plasma of human [14], complete loss of cleavage by a host furin [15], reduced binding of HLA to NF9 epitope, and reduced viral infectivity [16]. It has also been proven in cell culture and animal model studies that the mutation at the 614th position of the Spike causes high production of the pseudo-type virus particles and results in an increment of viral load [17].

In the presented manuscript, we have thoroughly identified and explored the mutability of the most important sites of the SARS-CoV-2 virus; i.e., the interface regions (residues 417–505) of Spike protein [18], as well as the relative influence of corresponding mutations on strengthening the interaction of the virus with host receptor in addition to detailing the mutational impact on the SARS-CoV-2 structural proteins and particularly on the Spike glycoprotein. To achieve the aforementioned goal we have used the sets of standard mutation modeling tools for identifying the Spike interface residues that, upon mutation, modify the affinity of ACE2 to Spike, and also the energy favorability of occurrences of these mutations. Further, MD simulation methods were applied to mutated Spike and interacting ACE2 proteins to find the net stability order of the Spike-ACE2 complex in consequences of Spike interface mutation. The pipeline of investigation applied in this research is presented in the flow diagram in Fig. 1 a, and the detailed process of mutational identification in Fig. S1.

Fig. 1

(a) Flow diagram representing computational methods and processes utilized to implement the current study. (b) Full-length trimeric Spike protein (individual chain colored in brown, light green, and yellow) interacting with ACE2 shown in light brown color. (c) Region of single-chain Spike interacting with ACE2 (encircled in black color) used in the mutational investigation. (d) Interface residues of Spike involved in interaction with ACE2 are shown in yellow color sticks conformation, while the red cartoon representation indicates the rest of the chain and the light brown cartoon indicates ACE2. (e) Interface showing positions of the WT residues of the Spike, which after getting mutated resides firmly at the interface (shown in orange), flees the interface (green), or could either reside or flee interface depending on a mutation (pale cyan). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Material and methods

SARS-CoV-2 Spike receptor-binding domain (RBD) complexed with human ACE2 was obtained from the protein data bank (PDB) with PDB ID: 6M17 [19], and the necessary protein preparation, including removal of small molecules/ligands, were done using PyMol software. The prepared protein complex was submitted to Protein Interface In Silico Mutation Scanning (PIIMS) tools for the evaluation of interface residues [20]. After evaluating the interface residues of the virus Spike protein, full length Spike-ACE2 complex PDB ID: 7DF4 [21] was prepared in PyMol and submitted to the mutation modeling tool Mutabind2 [22] for the evaluation of the structural impact of the mutations on the full-length Spike. For this purpose, each interface residue was mutated with the 19 other standard amino acids (AAs) among the known 20 proteinaceous AAs. Upon confirming the interface residues and associated mutations identified by Mutabind2, the protein complex was submitted to Dynamut2 for assessing the associated energetics, stability, and flexibility of the missense mutation [23]. Furthermore, to understand the impact of the mutations on the protein dynamics and relative stability strength, molecular dynamic (MD) simulation was performed using CABS-flex 2.0 [24]. The radius of gyration (ROG) of wild type (WT) and mutant proteins was also calculated using the supercomputer facility for Bioinformatics & Computational Biology, IIT Delhi [25]. At last, the associated secondary structure change position was evaluated in the Spike protein for the stabilizing mutation via the Multalin tool [26].

Results and discussion

Analysis of viruses for the identification of structural mutations at the AA level is highly beneficial not only in terms of cost-cutting but in identifying as early as possible the harmful variants that arise as a consequence of structural change. This is because of most of the genomic mutations observed in circulating SARS-CoV-2 viruses are expected to be either mildly harmful or neutral due to their encoding of neutral amino acids [11]. Therefore, to head-start the mutational investigation at the AA level in our process, Spike's interface residues were first identified using the PIIMS (Protein Interface In silico Mutation Scanning) server [20]. PIIMS functions in two modes: the first is to analyze the protein interface residues and hotspots of mutation, while the second is to analyze the effect of hotspot mutation.

Using the first mode, i.e. the Computational Alanine Scanning (CAS) method, the interface region of Spike protein, whose diagrammatic representation is shown in Fig. 1b and c, was targeted. CAS calculates interface residues based on the change in solvent-accessible surface area (ΔSASA) which is the difference of the SASA of protein subunit to their complex [27]. The other parameters obtained using CAS, in addition to interface residues, are binding energy and affinity fold change (shown in Table S1) which are the difference between wild type and mutational variant residue energy and affinity, respectively. The strength of mutation occurrence is related to the change in the binding energy (ΔΔGbind) [28], while the fold of affinity value indicates the hot-spot rank of the residue for being mutated.

The 21 Spike interface residues that were obtained via PIIMS CAS and were used in further mutational investigation are 417LYS, 445VAL, 449TYR, 453TYR, 455LEU, 456PHE, 473TYR, 484GLU, 486PHE, 487ASN, 489TYR, 490PHE, 493GLN, 494SER, 495TYR, 497PHE, 498GLN, 500THR, 501ASN, 503VAL, and 505TYR. These residues are labeled in yellow and in stick conformation in Fig. 1d, where the Spike protein chain is in red and the ACE2 receptor shown in light brown, while CAS-computed supplementary properties of these residues that are not required in our method of processing are shown in the Supplementary Table S1. After obtaining the details of Spike interface residues such as type and position by PIIMS, each residue position was replaced with 19 other standard amino acids sequentially using the Mutabind2 mutation modeling tools to observe the effect of the mutation on the interface of interaction. Therefore, a total of 399 mutations at 21 positions in Table S2 evaluated using Mutabind2 provided the output in terms of the binding affinity change (ΔΔGbind in kcal mol−1) with the details of mutated residue updated location with respect to the interface, i.e. its presence or absence on the interface and severity of mutation (deleteriousness).

Mutabind2 results help in scrutinizing the mutation presence/absence on the interface, and out of 399 evaluated mutations, 348 mutations at 20 interface positions were found to be sustained after getting mutated. Of these 20 positions, the 15 that were determined to be residing at the interface after getting mutation (enlisted here with WT residues) are 445VAL, 456PHE, 473TYR, 486PHE, 487ASN, 489TYR, 490PHE, 493GLN, 494SER, 495TYR, 498GLN, 500THR, 501ASN, 503VAL, and 505TYR, shown in orange sphere conformation in Fig. 1e. Meanwhile, 5 positions at the interface were identified as a gray region where certain mutations cause their disappearance from the interface, while the remainder causes their firm presence at the interface. These 5 positions given with WT residues are 417LYS, 449TYR, 453TYR, 455LEU, and 484GLU, shown in pale cyan sphere conformation in Fig. 1e.

At one position, namely position 497 (WT residue PHE), all 19 mutational variants fleeing their presence at the interface is shown in the green sphere in Fig. 1e. After identifying all interface mutations, the analytic tool Dynamut2 was employed to screen the behavior of these mutations, such as stabilization/destabilization and change in flexibility. Dynamut2 combines two modes: the first one uses a graph-based signature to represent wild type environment for investigating the effect of point mutation on protein stability and dynamics, while the second uses normal mode analysis to capture protein motion. The output of Dynamut2 provides the difference in the folding free energy of the WT protein to mutant variants; the positive values indicate a stabilizing nature of mutation while the negative value is indicative of destabilization. The stabilization/destabilization energy value of all the 348 interface mutations, of which 52 were found to be of stabilizing (+ΔΔGstability) and 296 were destabilizing free folding energy (-ΔΔGstability), are summarized in Table S3. The stabilization of free folding energy indicates the stability of the protein and it is one of the most important parameters to consider while investigating the mutations. Recently, researchers have also proven the dominant Spike mutation D614G running in the COVID-19 pandemic has a high Spike stability [29].

In the bar diagram of Fig. 2 a mutations of SARS-CoV-2 Spike interface residues are arranged as per their stabilizing energy, since the energy favourability of mutation occurrences is proportional to their change in free folding energy (±ΔΔGstability). In the given arrangement, the most highly favorable mutation obtained was Q498 M, while the least favorable was V503S. Further Dynamut2-assisted estimation of the changes of molecular interactions such as Hydrophobic interaction, Hydrogen bond, Ionic, Polar, and Vander wall interactions in WT and mutational variants determined that the reduced intramolecular interactions between the residues of SARS-CoV-2 indicate corresponding mutations which break the native interaction between residues, either to build the intermolecular interactions directly with ACE2 or to provide the necessary flexibility to the Spike protein for building the new interaction with ACE2.

Fig. 2

Effect of mutations on Spike stability; (a) Arrangement of Spike mutations in order of the free folding energy change (ΔΔGstability); (b) Stability derived from the RMSF values generated during the MD simulations. (c) Effect of Spike mutations on stability of ACE2. The stability order is derived from the RMSF of all ACE2 residues in MD simulation. (d) Effect of point mutations in Spike on the net stability of the Spike-ACE2. The net stability order is derived from the RMSF of all residues of the mutant Spike and RMSF of interacting ACE2 residues.

For instance, in Supplementary Table S4 comparing the mutations at residue position 417, the number of molecular interactions was reduced from 14 in WT to 10 for K417 M, to 11 for K417I and K417Y. In K417I and K417Y, where molecular interactions are the same in number, protein flexibility will be deciding factor for the preference of intermolecular bonding between Spike and ACE2. The alternation in the molecular interactions such as number of hydrophobic interaction, Hydrogen bond, van der Waals interaction, and Polar interaction as a consequence of point mutations at residue position 417 has been shown in Fig. 3 a to d, in which change in interaction forces was observed while comparing to the mutant variants. Similarly, the transition in the molecular interactions of all identified 52 energetically favorable Spike interface mutational variants is systematically displayed in Table S4.

Fig. 3

(a) Molecular interactions of the WT Spike protein residue Lysine 417 (shown in blue color). (b to d) Changes in the molecular interactions of the residue 417 (shown in the red color) of Spike caused by mutation to Isoleucine, Methionine, and Tyrosine, respectively. (e) RMSF obtained from the all-residue MD simulations of full-length WT and V503D mutated Spike (shown in black and red color respectively). (f) Zoomed-in RMSF profiles showing only interface residues of the WT and Mutant Spike. (g) RMSF obtained from the all-residue MD simulations of ACE2 complexed with the WT or V503D mutated Spike (shown in black and green, respectively). (h) Zoomed-in RMSF profiles showing only interface residues of ACE2 complexed to the WT and Mutant Spike. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Although the energy favorability value of mutations indicates the high chances of mutation occurrences, it doesn't reveal the severity of the mutation. Therefore, to investigate the mutational impact on interaction strength of Spike to ACE2, Molecular dynamic (MD) simulations have been performed using the CABSFLEX2.0 tool which is a coarse-grained protein modeling tool (CABS) based server for the fast and efficient simulation of protein flexibility. After obtaining the simulation results in the form of residues position versus RMSF of all residues of submitted structures containing interacting ACE2 (Chain A) to WT and mutated Spike (Chain B), the difference of the sum of all residues RMSF of chain A and B of mutated structure to the WT chain A and B were calculated using Equation (1).

Further percentage stability/destability was calculated by dividing the ΔRMSFx100 by the RMSF of WT using Equation (2). The calculated stability/destability of Spike in descending order as a consequence of mutations is shown in Fig. 2b, indicating the highest stability of Spike in isolation was conferred by mutation V503D, while the lowest stability was caused by V503 N. In Fig. 2c the change stability orders of interacting ACE2 in consequences to Spike interface mutation are shown too, the given figure indicates the highest stability of ACE2 is convened by mutation Q493 M, while the least is by V503H. The net stability of mutated Spike-ACE2 complexes was calculated by the addition of the stability of each component, i.e. percentage stability of Spike and ACE2 stability. The obtained results are arranged in descending order and shown in Fig. 2d indicating the net stability of the complex of Spike-ACE2 as a consequence of different point mutations, in which 36 lead to enhanced stability of Spike-ACE2 while 16 destabilize the complex.

The 36-point mutations that increase the stability of the Spike-ACE2 complex, and the percent increase in the stability in comparison to the WT complex, are shown in Table 1 . The flow diagram Fig. S1 explains the identification process of 36 mutational variants that highly stabilize the Spike-ACE2 and could be potentially harmful, in the very simplest manner. Among the 36 stabilizing mutations, the highest stabilization mutation obtained was V503D, which is also known to be neutralization-resistant [30], while the second-ranked mutation Q493 M is proven to significantly enhance the Spike-ACE2 interaction [31,32]. The third-ranked mutation in our analysis E484 M is known for spreading faster in Brazilian variants [33]. Identified Spike mutations in our analysis K417Y and E484A were recently reported to be a variant of concern for SARS-CoV-2:Omicron that requires the attention of a vaccine developer [34]. In Table S5, all identified 36 novel mutational variants that could be potentially harmful are shown. The published outcomes for the impact of a few of these identified variants by different researchers are also referenced in the given Table S5.

Table-1

Energetically favorable (MD simulation-based) stabilizing mutation of Spike-ACE2 complex, their ROG, and percentage increase in the stability of mutated Spike-ACE2 complex.

S.No.	Mutation	ROG IN Å	% Increase in Stability of Spike-ACE2 complex
1	V503D	25.16	29.4
2	Q493 M	25.16	26.8
3	E484 M	25.16	22.8
4	K417I	25.16	19.8
5	F490C	25.16	19.0
6	Q498Y	25.15	17.6
7	Y505R	25.16	17.2
8	S494 M	25.16	16.9
9	K417Y	25.16	15.9
10	E484W	25.15	14.5
11	Q498I	25.16	13.3
12	V503S	25.16	12.8
13	S494P	25.16	11.7
14	F486Y	25.16	10.7
15	Y449K	25.15	10.6
16	E484H	25.16	10.0
17	Y505K	25.16	9.8
18	V503H	25.16	9.5
19	V445D	25.16	8.1
20	V503G	25.16	6.8
21	V503E	25.16	5.6
22	Q498 M	25.16	5.1
23	Q498L	25.16	4.3
24	E484R	25.15	3.6
25	Y449H	25.15	3.5
26	Q493I	25.16	3.5
27	S494V	25.16	2.6
28	S494K	25.16	2.5
29	E484A	25.16	2.4
30	Y449R	25.15	2.3
31	S494I	25.16	1.4
32	E484L	25.15	1.3
33	V503P	25.16	1.2
34	S494L	25.16	1.1
35	Q498V	25.16	0.5
36	Q493L	25.16	0.1

Full-length Spike and interacting ACE2 residues' root mean square fraction (RMSF) of highly stabilizing mutational variants are shown in Fig. 3e and g, and although noticeable differences can be visualized the deviation differences are not very clear. Therefore, Fig. 3f and h shows reduced residues’ (only interacting residues) RMSF, confirming the stability differences between WT and mutant proteins. Further, to confirm whether the stabilization of mutated Spike protein is caused by shrinkage of its 3D structure or the altered interactions occur as a result of the position of secondary structure change (SSC), the calculation of Radius of Gyration (ROG) and SSC analysis of all 36 mutations that stabilize the Spike-ACE2 complex were performed, and a very minimal change in ROG in Table 1 and SS in Table S5 were obtained, indicating that mutated protein stability (or reduced fluctuation in RMSF) is not affected by ROG. The changed secondary structure (SS) locations revealed the importance of these positions for the stability of the complex of mutated Spike to ACE2. Our analysis of SS change as an impact of Spike point mutations reveals that very little change in SS is in good agreement with previous studies [35,36].

Conclusions

In conclusion, 36 novel mutations that highly stabilize the Spike-ACE2 interaction were identified. First, the interface position and associated residues of Spike protein were identified, and further mutational variants that reside on the interface were discovered. The Spike-ACE2 destabilizing mutations were filtered from the stabilizing ones using the mutational modeling tools. The detailed interaction profile and energy favorability of Spike stabilizing mutations are explored at this stage. At the next stage, MD simulation for Spike stabilizing mutations’ was performed to determine the relative stability order for the Spike-ACE2 complex. The analysis of the order of the net stabilities of the Spike-ACE2 complexes as the consequences of the point mutations in Spike revealed that 36 mutations are expected to stabilize the complex, while 16 mutations can destabilize it. The presented mutational impact on the strengthening of Spike-ACE2 interaction and their feature of change in the molecular interaction and secondary structure point towards the preferred ordered basis of exploring the mutation-driven unexplored mechanism in viral pathology. In addition to providing the sequential basis for mechanistic exploration, the presented outcomes will help in the development of potent vaccines and novel antivirals against SARS-CoV-2 variants, in addition to finding the solution for counteracting the future consequences of the SARS-CoV-2 outbreak.

Declaration of competing interest

The authors declare no conflict of interests.