SARS-CoV-2 spike evolutionary behaviors; simulation of N501Y mutation outcomes in terms of immunogenicity and structural characteristic.

Since the emergence of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), a large number of mutations in its genome have been reported. Some of the mutations occur in noncoding regions without affecting the pathobiology of the virus, while mutations in coding regions are significant. One of the regions where a mutation can occur, affecting the function of the virus is at the receptor-binding domain (RBD) of the spike protein. RBD interacts with angiotensin-converting enzyme 2 (ACE2) and facilitates the entry of the virus into the host cells. There is a lot of focus on RBD mutations, especially the displacement of N501Y which is observed in the UK/Kent, South Africa, and Brazilian lineages of SARS-CoV-2. Our group utilizes computational biology approaches such as immunoinformatics, protein-protein interaction analysis, molecular dynamics, free energy computation, and tertiary structure analysis to disclose the consequences of N501Y mutation at the molecular level. Surprisingly, we discovered that this mutation reduces the immunogenicity of the spike protein; also, displacement of Asn with Tyr reduces protein compactness and significantly increases the stability of the spike protein and its affinity to ACE2. Moreover, following the N501Y mutation secondary structure and folding of the spike protein changed dramatically.

INTRODUCTION

The emergence of COVID‐19 can be considered as one of the main human health challenges in the present century. COVID‐19 is caused by severe acute respiratory syndrome coronavirus 2 (SARS‑CoV‑2), which is a positive‐strand RNA virus. The genome length of this virus is 29 903 nucleotides (Wuhan strain) and capable of encoding 12 structural and 16 nonstructural proteins. In addition to the high rate of transmission, mutations in genome‐encoding regions have made it difficult to control the SARS‑CoV‑2 pandemic. Mutation in all SARS‑CoV‑2 proteins can be significant; however, amino acid displacement in spike protein is important due to its role in binding and the virus entering into the host cells through interaction with the angiotensin‐converting enzyme 2 (ACE2). Although not much time has passed since the onset of COVID‐19, many mutations and polymorphisms have been reported in the gene encoding the spike protein. Mutations and polymorphisms as an unfortunate consequence of evolutionary pressures, , , , based on their position in the spike sequence they are divided into two categories; inside the ACE2 binding interface and outside this region. Both types of mutations can affect the stability and affinity of the spike to its receptor. An important mutation in the ACE2‐binding domain is N501Y, which has been observed in the UK/Kent (B.1.1.7; N501Y.V1), South Africa (B.1.351; N501Y.V2), and Brazilian (B.1.1.28; P1; N501Y.V3) spike protein variants, generating the explanation of variants of concern and suggestion for lineage‐specific surveillance.

So far, many studies have evaluated the biological consequences of the N501Y mutation in SARS‐CoV‐2 behavior. In an investigation which 496 patients were enrolled to evaluate the difference between clinical signs and genetic characteristics of B.1.1.7 versus non‐B.1.1.7 variants, it was found that viral load was higher in B.1.1.7. These samples have a lower cycle threshold in a quantitative polymerase chain reaction assay, and the read depth was higher in B.1.1.7 carrying samples. Examination of other variants carrying the N501Y mutation, such as B.1.351 and B.1.1.28, have shown that the ability to escape from neutralizing antibodies responses increases in these lineages. , In addition to characterized lineages, the prevalence of unknown mutants of SARS‐CoV‐2 carrying the N501Y has been reported and raised concerns. For example, a clinical study that was conducted to investigate the prevalence of COVID‐19 in Switzerland identified many patients infected by incognito mutants of spike that did not belong to any lineage. SARS‐CoV‐2 linages that introduced K417N/E484K/N501Y spike mutations show a greater transmission and infectivity potency compared to the native form. The ability of the Pfizer vaccine to immunize against variants containing these three mutations is reduced by about 6–8‐fold. Liu et al. in an investigation that evaluated the transmission capability of the B.1.1.7 variant, they realized that N501Y substitution is the main factor in increasing the rate of transmission and spike affinity to its receptor. Also, spikeN501Y exhibits a high ability to neutralize the therapeutic effects of antibodies composed of V‐region IGHV3‐53. Regarding the importance of the N501Y mutation in the spike protein biological behavior, our group studied the structural characteristics of this variant using computational biology and structural bioinformatics approaches, such as immunoinformatics, docking, molecular dynamics (MD), and free energy computations. Since limited numbers of studies were conducted to analyze the three‐dimensional structure of spikeN501Y, the results of this study can provide an atomic‐level information for biology researchers and other related disciplines.

MATERIALS AND METHODS

Immunogenicity analysis

The N501Y mutation is located in the interaction interface of spike protein with its receptor. Clinical studies demonstrate that this replacement can affect the neutralization resistance potency of SARS‑CoV‑2 and its viral load, which may indicate a change in immune responses against this variant. Immunogenicity of the N501Y variant was compared with native spike using the immune epitope database (IEDB) server which is funded by the National Institute of Allergy and Infectious Diseases. , This server is a comprehensive immunoinformatics database for examining T‐cells and B‐cells antigens. The sequences of native and N501Y spike variants were retrieved from UniProt in FASTA format. Then, for identification and scoring of T‐cells epitopes, major histocompatibility complex Class I (MHC‐I) and MHC‐II binding predictions of IEDB were implemented. For this purpose, IEDB recommended 2.22 (MHC‐I) and NetMHCPan 4.1 EL (MHC‐II) methods were chosen. As well Linear and discontinuous B‐cells epitopes were predicted by Antigen Sequence Properties (Bepipred Linear Epitope Prediction 2.0) and DiscoTope tools of IEDB, respectively. ,

Docking and residue scanning methodology

Molecular docking is a versatile and high throughput approach for molecular interactions modeling. Docking and related items such as virtual screening have revolutionized the world of drug design and biomolecule studies. Due to the importance of the subject, many algorithms and docking tools have been developed. The first category is global algorithms (also called blind docking) which docking are performed without accurate information about the binding pockets, like the implemented algorithm in the MDockPeP server. Another category is called local docking algorithms and requires sufficient information from the binding pocket, like the algorithm employed on the HADDOCK server. Some tools also provide a hybrid environment that allows performing both types of docking such as Rosetta, Schrodinger, ClusPro, and AutoDock. In the present study, since the residues in the spike–ACE2 interface was previously identified, local docking was performed by Rosetta v3.12. Throughout docking with Rosetta, parameters such as rotation, translation, and the number of outputs (‐nstruct) were set at 8°, 3 Å, and 1000, respectively. In the present study, desired structures of native (Protein Data Bank [PDB] ID: 7DF4) and mutant (PDB ID: 7MJN) spike–ACE2 complexes , were retrieved from RCSB (http://www.rcsb.org). All input structures before docking conducted to the relaxation and repack process for backbone and sidechains clashes fixation, respectively. Furthermore, a local refinement protocol was done for local docking outputs to improve predicted quality poses. In parallel with docking, residual scanning analysis was performed using DynaMut server to investigate the effect of the mutation at position 501 on the stability and affinity of the spike for its receptor. The highest score and best conformations of the spike–ACE2 complex obtained in docking were conducted to other analysis steps.

MD simulation

Various laboratory techniques such as X‐ray, nuclear magnetic resonance, and cryo‐electron microscopy examine changes in the secondary and tertiary structure of proteins. , Despite the high accuracy of these approaches, some disadvantages, such as complexity and high cost limit their use. In parallel with these techniques, computational biology has advanced significantly in recent years. Affordability, high speed, and flexibility are the main advantages of computational and in‐silico methods. One of the most popular techniques in the computational biology of biomolecules is MD simulation (also called computational molecular microscope) which enables investigators to determine changes at the atomic level. In MD, the behavior of molecules over time is investigated using Newtonian equations. Our group utilizes MD to investigate the effect of N501Y displacement on the structural characteristics of the spike. All MD simulation steps were performed using GROMACS v2021 (CUDA compiled) on Linux operating system (Ubuntu v20.04 LTS). CHARMM36 all‐atom force field was implemented to parameterize the structures and generation of topologies. To equilibrate systems, the Canonical ensemble (constant number, volume, temperature) was performed for temperature coupling at 310 K by V‐rescale algorithm for 150 ps, as well as isothermal–isobaric ensemble (constant number, pressure, temperature) by Parrinello–Rahman algorithm for 500 ps was accomplished for pressure coupling at 1 bar. , During equilibration, constraints were applied for h‐bonds; furthermore, rcoulomb and rvdw were set to 1.2 nm. The TIP3P water model was used to solvation of systems and Na+ and Cl­ ions at a concentration of 0.15 M were added for neutralization. In the MD production step, all systems were simulated for 30 ns. Finally, structural terms such as root‐mean‐square deviation (RMSD), root‐mean‐square fluctuation (RMSF), solvent accessible surface area (SASA), the radius of gyration (Rg), free energy landscape (FEL), principal component analysis (PCA), secondary structure, hydrogen bonds, minimum distance, and the number of contacts were extracted from the output trajectory files.

Free energy computation

Investigating the affinity and stability, following mutations is a fundamental step in the study and design of proteins. One of the main ways to evaluate these biophysical parameters is free energy calculation. Although numerous algorithms including free energy perturbation and thermodynamic integration were developed for this purpose, low speed limits their use to small molecules studies. Another available method for calculating free energy is molecular mechanics Poisson−Boltzmann surface area (MM‐PBSA), which has high speed and accuracy. This algorithm decomposes free energy terms into potential energy and solvation energy. The potential energy consists of both bonded and nonbonded interactions energies and is computed using MM functions. In addition, the solvation energy term is divided into polar solvation energy (electrostatic term) and nonpolar solvation energy (nonelectrostatic term). To calculate the binding free energy of the native and N501Y spike–ACE2 complexes, g_mmpbsa which is equipped with the MM‐PBSA algorithm was implemented. During free energy calculation, the SASA method was chosen for nonpolar solvation energy computation and ‐incl_14 set at yes. Besides, we employed Rosetta (scoring protocol), FoldX, and OSPREY (v3.2.2) for the structure extracted from the last frame of the simulation trajectory to computation structural energy of native and mutant spike in unbounded mode. OSPREY and FoldX are versatile tools for tracking mutation in proteins. ,

Tertiary structure topology analysis

To rigorously investigate the impact of the N501Y mutation on the folding characteristics of spike protein, structural alignment along with three functions including the probability density function (PDF), contact map, and tunnel analysis were applied for the residues around the mutated position. PDF was examined based on dihedral angle and triangle area using G_Measure (v.0.8); in addition, contact maps for carbon alpha and tunnel analysis were analyzed by utilizing Pymod (v3.0) and Caver (3.0.3), respectively. ,

RESULTS AND DISCUSSION

Insight into immunological characteristics of native and N501Y variants of spike comprehensively was obtained in two steps. First, MHC‐I and MHC‐II epitopes were predicted for the human leukocyte antigen allele reference set. It was observed that the N501Y spike has the less binding ability for MHC‐1 compared to the native type. Then, linear and discontinuous epitopes for B‐cells were anticipated. The mutation at position 501 significantly reduced the immunogenicity of adjacent residues which forms the main epitope for B‐cells (Table 1).

Table 1

IEDB server outputs, as seen native spike exhibited a higher immunogenicity potency for immune cells (higher score indicates higher affinity)

Position	DiscoTope	ASP	MHC‐I		MHC‐II
498W	1.715	0.574	N	0.484	N	69
498M	0.196	0.561
500W	3.006	0.584
500M	1.737	0.558
501W	−0.797	0.579
501M	−1.935	0.550	M	0.483	M	69
502W	1.204	0.585
502M	0.107	0.560
505W	−2.396	0.563
505M	−3.248	0.535

Abbreviations: ASP, antigen sequence properties; IEDB, immune epitope database; M, mutant; MHC‐I, major histocompatibility complex Class I; N, native.

In addition to the higher binding affinity, the number of epitope regions in the native spike was higher than the N501Y variant (Figure 1). In the native form of a spike, 23 residues were identified as epitopes, while the value for the mutant variant was 21 residues. In previous studies, , , confirmation of lineages carrying the N501Y have higher viral load and resistance to antibody neutralization, the observed decrease in immunogenicity can be attributed to the occurrence of this phenomenon.

Figure 1

Structural epitopes examined using DiscoTope (same and different epitopes are shown in cyan and hot pink, respectively)

Docking and residue scanning

The pattern of residues interactions for spike–ACE2 complex in both native and N501Y variants was tracked using Rosetta docking protocol as a state‐of‐the‐art, free, and open‐source suite. Rosetta docking energy calculates using physics–statistics functions and is reported in the form of a Rosetta energy unit. For native spike, Asn at position 501 forms a hydrogen bond (2.71 Å) with Lys353. Furthermore, Thr500 interacts with residue Asp355 with two hydrogen bonds (2.61 and 3.06 Å) versus in the N501Y variant, Tyr at position 501 is able to form one hydrogen bond with Lys353 (2.93 Å), while its aromatic ring forms two hydrophobic bonds with Tyr41 (5.13 Å) and Lys353 (3.94 Å) ACE2. In addition, Thr500, like the native spike, forms two hydrogen bonds with Asp355, but with a shorter distance (both are 2.63 Å). Rosetta energies for native and N501Y variants were −1945.63 and −1948.516, respectively. Indeed, following the N501Y mutation affinity of spike protein increases for its receptor (Table 2).

Table 2

Some energy terms reported by Rosetta (the unit of all energy terms is REU)

Spike	fa_atr	fa_elec	hbond_sr_bb	hbond_lr_bb	hbond_bb_sc	hbond_sc	Omega	Total
Native	−4950.7	−1242.1	−307.867	−87.296	−138.802	−104.010	82.678	−1945.57
Mutant	−4949.6	−1242.2	−307.503	−88.195	−141.318	−101.883	83.503	−1948.51

In addition to the distance and number of interactions between the residues in the spike interface with its receptor, another factor that increases the binding potency in the N501Y is the shape of the Tyr side chain. Unlike Asn501 which is adjacent to Lys353, the side chain and aromatic ring of the Tyr501 are surrounded by Y41 and Lys353 (Figure 2). This characteristic leads to the better formation of the binding pocket of the N501Y SARS‐CoV‐2 lineage. Residue scanning calculations of DynaMut also confirm the docking results and show that the stability of the spike–ACE2 complex improves following the N501Y mutation (−0.522 kcal/mol).

Figure 2

Rosetta docking outputs for native and N501Y spike variants

MD investigations

Investigation of system equilibration in terms of thermodynamic parameters is an essential step before analyzing the MD outputs. For this purpose, some thermodynamics parameters such as pressure, temperature, and density were extracted for the studied systems. Examination of the mentioned parameters indicated that the target systems were equilibrated (Table 3).

Table 3

As can be seen, all variables coupled to desired conditions

System	Temperature (K)	Pressure (bar)	Density (kg/m3)
Native	309.99	1.2	1044.4
Mutant	310	0.99	1044

To accurately assess the equilibration and to make sure that MD time was sufficient, RMSD curves were extracted for the target systems. We found that during the simulation, this structural parameter also reached equilibrium (Figure 3).

Figure 3

Root‐mean‐square deviation (RMSD) curve for the native and mutant spike, after 18 ns, the amount of conformational shift in both systems has been minimized

To investigate the role of amino acids in the rate of structural changes, the RMSF factor that indicates residue fluctuations was assigned. As can be seen in Figure 4, the number of structural fluctuations following the mutation at position 501 decreases, this is consistent with the results of RMSD output. Replacing Asn with Tyr reduces deviations and variations in the RBD domain of spike. The occurrence of this phenomenon can be related to the interaction potency of Tyr. Since Tyr compared to Asn has a higher ability to interact with adjacent residues in the RBD domain; thus, it applies less fluctuation to the spike structure and can improve the stability of this fragment of the spike. The mean fluctuations in the native and mutant structures were 0.28 and 0.26 nm, respectively.

Figure 4

 Root‐mean‐square fluctuation (RMSF) diagram, the higher number of polar interactions in the mutant compared to the native spike (4 vs. 2) reduces the structural fluctuations rate following the N501Y mutation

The effect of the N501Y mutation on spike folding was tracked by evaluating the Rg and SASA functions. At first, the Rg that referred to the protein compactness was measured in both native and mutant species. Rg represents the average distance of the atoms that make up a molecule from its center. The average Rg was 0.53 nm in the native and 0.54 nm in the mutant. This finding indicates a reduction in the compactness of the spike protein following mutation. In addition, the average of SASA for native and mutant spike were 11.4 and 11.7 nm2, respectively (Figure 5). The results of the Rg and SASA investigations were highly consistent. This signifies that placement Tyr at position 501 has a great impact on the folding of adjacent areas. The increase in accessibility and interaction interface observed in mutant variants can be one of the main factors in promoting spike affinity to ACE2 receptors.

Figure 5

(A) Radius of gyration (Rg) and (B) solvent accessible surface area diagram for both native and mutant spike

In the next step, FEL and PCA functions were extracted to track the effect of the N501Y mutation on the folding and the conformational energy patterns of the target structures (Figure 6). Following the mutation, spike folding significantly affects the conformational energy patterns in the mutant species, where it was more concentrated compared to the native form. These results confirm the RMSD, RMSF, and Rg outputs and indicate that the flexibility and structural variation decrease following N501Y mutation.

Figure 6

Free energy landscape and principal component analysis (PCA) diagram. Rg, radius of gyration; RMSD, root‐mean‐square deviation

Analysis of secondary structure elements using the defined secondary structure of proteins module of GROMACS (residues 490–510) was another examined parameter. We observed that the N501Y displacement could affect the composition of the secondary structural components of its adjacent amino acids. Compared to the native structure, in the mutated spike amount of bend and 3‐helix were decreased, while the coil elements frequency had increased (Figure 7).

Figure 7

Schematic representation of secondary structure for the native and mutant spike

In the final step, MD trajectory files were examined to track the effect of N501Y mutation on the affinity of spike for its receptor using functions such as a number of hydrogen bonds, the distance between interface residues, and the number of contacts at distance less than 0.6 nm for the spike–ACE2 complexes. Hydrogen bond analysis revealed that minimum, maximum, and the average of hydrogen bonds for complex of native spike with ACE2 were 0, 6, and 2.95, respectively. However, these values for mutant spike were 1, 8, and 4.7 (Figure 8).

Figure 8

Hydrogen bonds state throughout molecular dynamics simulation

Evaluation of distance between Asn and Tyr at position 501 with interface residues of ACE2 showed that the mean distance in mutant species (0.111 nm) was less than the native (0.112 nm), which could be evidence for greater ability of spike to ACE2 cell surface receptor. Finally, contacts at distance less than 0.6 nm for native and mutant complexes were investigated and the results were significantly consistent with other observations. We found that the number of contacts for mutant complex (886) is higher than native (755) lineage (Figure 9). These results along with the output of hydrogen bonds investigation confirm the higher binding affinity in mutant spike for ACE2 compared to native form.

Figure 9

(A) Number of contacts and (B) minimum distance for interface residues in spike–angiotensin‐converting enzyme 2 complexes

Regarding the results observed during MD simulation, it can be claimed that the N501Y mutation has a significant effect on the RBD domain of SARS‑CoV‑2 spike protein and can increase the spike stability and its affinity to ACE2.

Free energy analysis

Energy calculation was implemented for two aims. First, the affinity of native and mutant spike for ACE2 was investigated using binding free energy approximation. The findings of this step confirmed the results of docking and MD. We observed that the binding affinity for the ACE2 increases following the replacement of Asn with Tyr (Figure 10).

Figure 10

Schematic representation of the binding energy for native and N501Y carrying complexes

Further examination was done by decomposition of binding free energy. We found that, with the exception of polar solvent energy, the other constituent functions of the binding energy were more optimal for the complex containing the N501Y mutation (Table 4).

Table 4

Binding free energy decomposition for spike–ACE2 complexes

Complex	VdW energy	ES energy	PS energy	SASA energy	Binding energy
Native	−283.415	−2353.393	533.526	−37.414	−2140.696
Mutant	−313.698	−2487.918	593.290	−39.733	−2248.058

Abbreviations: ACE2, angiotensin‐converting enzyme 2; ES, electrostatic energy; PS, polar solvation energy; SASA, solvent accessible surface area; VdW, Vander walls energy.

Furthermore, per residue binding free energy decomposition was carried out for residues 41, 353, and 355 ACE2 and residues 500–502 spike. As expected, compared to Asn501, placement of Tyr imposes much higher binding energy to the spike for binding to ACE2 (Figure 11).

Figure 11

Binding free energy (kJ/mol) of position 501 and surrounding amino acids for spike and their interface residues in angiotensin‐converting enzyme 2

After studying the binding free energy, to accurately examine the effect of the mutation at position 501 on the stability of the spike protein, we investigated the folding energy of this protein in the unbounded form. For this purpose, FoldX, Rosetta, and OSPREY were employed. Examination of the results of FoldX, Rosetta, and OSPREY showed that the N501Y mutation increases the stability of the spike protein (Table 5). The outputs of the structural energy investigation agreed with the results of other analysis steps and indicate the positive effect of the N501Y on the stability and affinity of the spike to its receptor.

Table 5

Stability assessment of spike in native and mutant forms in terms of folding energy

Lineage	Tools
	FoldX(kJ/mol)	Rosetta(REU)	OSPREY(kJ/mol)
Native	−117.78	−382.720	−24.76
Mutant	−141.9	−402.669	−38.56

Abbreviation: REU, Rosetta energy unit.

During structural alignment, we observed a structural shift (2.2 nm) in the structure of the mutant spike compared to the native structure, which is a reason for the change in the binding affinity of the spike to interaction with ACE2 (Figure 12).

Figure 12

Structural alignment of native (red) and mutant (cyan) spike

PDF analysis consistent with results of MD and structural alignments revealed an extensive change in the folding of the spike. We found that the pattern of dihedral angles and triangles area changes following the N501Y displacement. Furthermore, analysis of the contact map for carbon alpha at position 501 showed that the arrangement of residues around this position undergoes extensive changes, in the mutant structure Tyr501 interacts with residues 497–506, while in the native spike, Asn501 interacts with residues 498–506 (excluding position 504). More stable interaction of position 501 with surrounding residues in mutant spike maybe reduce the structural fluctuations and increase its stability. Besides, the study of the tunnels around the mutated position revealed there is more space in the mutant spike that can facilitate the interaction of spike with ACE2 (Figure 13). The Caver results showed that the length and radius of the available tunnel around the Tyr501 (2 and 1.2 Å) are greater than the Asn501 (1.9 and 1.1 Å).

Figure 13

Representation of a probability density function, contact map, and tunnel analysis outputs for desired structures

CONCLUSION

One of the main challenges in detecting and controlling SARS‐CoV‐2 is mutation occurring in functional regions of the virus genome, such as the coding sequence of the spike protein. So far, many mutations (E484K, Q493N, and N501T) have been identified in the RBD domain of this protein. Given that this protein is the main target in vaccination against COVID‐19, tracking the effect of these mutations, and changing the function of the spike is very valuable. Our team in this current study uses various computational biology approaches examining the significant effects on the function of the spike molecular behavior due to the N501Y mutation. We discovered that this mutation reduces immunogenicity, increases stability, and affinity for the spike receptor.

Since this mutation can facilitate the potency of the virus entering the cell and exacerbate its clinical symptoms. Due to the N501Y mutation effect on the immunogenicity of the virus induces immune responses, we recommend that more attention be paid to the identification of SARS‐CoV‐2 lineage during the admittance of COVID‐19 patients. Having efficient diagnostic and therapeutic protocols implemented for patients carrying this mutation compared to the reference lineage (Wuhan strain).

CONFLICT OF INTERESTS

The authors declare that there are no conflict of interests.

AUTHOR CONTRIBUTIONS

Methodology and software: Neda Rostami. Analyzed data and co‐wrote the paper: Edris Choupani. Software and co‐wrote the paper: Yaeren Hernandez. Performed bioinformatic analyses: Seyed S. Arab. Conceptualization: Seyed M. Jazayeri. Designed experiment, data analysis, and conceptualization: Mohammad M. Gomari.