N439K Variant in Spike Protein Alter the Infection Efficiency and Antigenicity of SARS-CoV-2 Based on Molecular Dynamics Simulation

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), causing an outbreak of coronavirus disease 2019 (COVID-19), has been undergoing various mutations. The analysis of the structural and energetic effects of mutations on protein-protein interactions between the receptor binding domain (RBD) of SARS-CoV-2 and angiotensin converting enzyme 2 (ACE2) or neutralizing monoclonal antibodies will be beneficial for epidemic surveillance, diagnosis, and optimization of neutralizing agents. According to the molecular dynamics simulation, a key mutation N439K in the SARS-CoV-2 RBD region created a new salt bridge with Glu329 of hACE2, which resulted in greater electrostatic complementarity, and created a weak salt bridge with Asp442 of RBD. Furthermore, the N439K-mutated RBD bound hACE2 with a higher affinity than wild-type, which may lead to more infectious. In addition, the N439K-mutated RBD was markedly resistant to the SARS-CoV-2 neutralizing antibody REGN10987, which may lead to the failure of neutralization. The results show consistent with the previous experimental conclusion and clarify the structural mechanism under affinity changes. Our methods will offer guidance on the assessment of the infection efficiency and antigenicity effect of continuing mutations in SARS-CoV-2.

Introduction

Coronavirus disease 2019 (COVID-19), caused by a single-stranded positive-strand RNA virus named severe acute respiratory syndrome-coronavirus 2 (SARS-CoV-2), is a major threat to public health worldwide (Wang et al., 2021). As of 6 June 2021, over 172 million confirmed cases including more than 3.7 million COVID-19 related deaths were reported worldwide according to the World Health Organization.¹ It has been reported that SARS-CoV-2 infects humans through the binding of the homo-trimeric spike (S) glycoprotein to human angiotensin converting enzyme 2 (hACE2), and this infection mechanism for viral entry is also used by SARS-CoV (Li et al., 2003; Zhou et al., 2020). The surface spike glycoprotein including S1 and S2 subunits is the major antigen of coronaviruses, and S1 binds to host cells whereas S2 mediates viral membrane fusion. The receptor-binding domain (RBD) mediates the binding of the virus to host cells, which is a critical step for viral entry (Hoffmann et al., 2020; Lu et al., 2020; Walls et al., 2020; Zhou et al., 2020). According to the high-resolution crystal structure (Lan et al., 2020; Yan et al., 2020), the receptor-binding motif (RBM) is essential for RBD and contacts highly with hACE2. The structural characterization of the pre-fusion S protein provides atomic level information to guide the design and development of antibody (Wrapp et al., 2020; Wu et al., 2020; Wang et al., 2021).

Neutralizing monoclonal antibodies of the immune system, which play an important role in fighting against viral infections, have been found to target the SARS-CoV-2 RBD and exert neutralization activity by disrupting the virus binding (Cao et al., 2020; Hansen et al., 2020; Ju et al., 2020; Li F. et al., 2020; Liu et al., 2020; Shi et al., 2020). During the virus transmission, alterations of amino acid in the surface spike protein may significantly alter the virus antigenicity and the efficacy of neutralizing antibodies. As SARS-CoV-2 spreads around the world, mutations in spike protein had been continuously reported (Korber et al., 2020; Morse et al., 2020; Saha et al., 2020; Sheikh et al., 2020; van Dorp et al., 2020). There are 930 naturally occurring missense mutations in SARS-CoV-2 spike protein that had been reported in the GISAID database (Supplementary Table 1), and a key mutation from ASP614 to GLY614 (D614G) in SARS-CoV-2 spike protein confer the SARS-CoV-2 more infectious than the original strain (Li Q. et al., 2020). However, most vaccines, testing reagents, and antibodies for SARS-CoV-2 are based on the S protein of the Wuhan reference strain (GenBank: MN908947.3). Among other coronaviruses, missense mutations had been demonstrated to confer resistance to neutralizing antibodies in MERS-CoV and SARS-CoV (ter Meulen et al., 2006; Sui et al., 2008; Tang et al., 2014). In the case of HIV, missense mutations are known to influence envelope glycoprotein expression, virion infectivity (Asmal et al., 2011), alter neutralization sensitivity (LaBranche et al., 2019), and confer resistance with neutralizing antibodies (Bricault et al., 2019; Zhou et al., 2019). The epidemiological observations have proved the mutations, especially in S protein, provide a plausible mechanism for the increased observed infectivity and bring challenges to antibody development for SARS-CoV-2.

Here, we evaluated genomic sequences of SARS-CoV-2 and observed a high frequency (0.72%) mutation occurring in the RBM region, N439K, which was first sampled in March 2020 in Scotland from lineage B.1 on the background of D614G, has arisen independently multiple times. It was observed in 34 countries and the second most commonly observed RBD mutation globally. Study have demonstrated that N439K S protein enhances binding affinity to the hACE2 receptor and reduces the neutralization activity of monoclonal antibodies (Thomson et al., 2021). To further provide insights into the infection and neutralization processes of N439K at the molecular level, we present performed molecular dynamics (MD) simulations of the binary complexes of the RBD domain with the common receptor hACE2 and the neutralizing monoclonal antibody (mAb) CB6/REGN10987, respectively. The structure of N439K-mutated RBD-hACE2 complexes show a new salt bridge and local interaction. The energetic details further indicate that the mutated RBD-hACE2 complexes show higher affinity than the original strain which are mainly attributed to the changes in van der Waals and electrostatic energies of some key residues. On the other hand, N439K reduced the sensitivity to neutralizing antibodies which are mainly attributed to polar solvation and electrostatic interactions. Taken together, the SARS-CoV-2 spike protein with N439K may be more infectious and become resistant to some SARS-CoV-2 neutralizing antibodies.

Materials and Methods

Multi-Sequence Alignment and Structure Preparation

In this study, SARS-CoV-2 sequences were aligned against the Wuhan reference genome (GenBank: MN908947.3) using MAFFT (Katoh and Standley, 2013). The wild-type crystal structures, including RBD-hACE2 complex (PDB:6M0J), RBD-CB6 complex(PDB: 7C01) and RBD-REGN10987 complex(PDB: 6XDG), were directly downloaded from the Protein Data Bank² and removed the water molecules using VMD (Humphrey et al., 1996) software. For mutated complexes, firstly wild-type hACE2 and mAbs structures were obtained by removing the RBD structures from RBD-ACE2 and RBD-mAbs complexes, respectively. Then the mutated SARS-CoV-2 RBD (N439K) structures were constructed using homology modeling by SWISS-MODEL (Waterhouse et al., 2018). The mutated SARS-CoV-2 RBD was aligned with hACE2 and mAbs based on the RBD-hACE2 and RBD-mAbs structures using PyMOL software. All structural figures were generated utilizing PyMOL software (Lilkova et al., 2015). Finally, the Ramachandran Plot (Laskowski et al., 1993), ANOLEA (Melo et al., 1997), Z-Score (Pontius et al., 1996), Verify 3D results (Bowie et al., 1991), PROCHECK (Laskowski et al., 1993) and Q-Mean (Waterhouse et al., 2018) method were used to check the quality of model.

MD Simulation

For each system, MD simulations were performed using GROMACS 5.1.4, the latest CHARMM36 (Tortorici and Veesler, 2019)force field was selected, an explicit solvent model was used by an explicit TIP3P model (Price and Brooks, 2004). Then complexes were solvated in a rectangular periodic box and had a 10 Å buffer distance from along each side, then a suitable number of Na+ and CL- ions were added to neutralize the whole system and mimicked a salt solution concentration of 0.15M. Specifically, the following stages of MD simulations were performed before the production simulation: (1) A full 50,000-step energy minimization was performed with all atoms unrestrained. (2) Each complex was equilibrated with NVT (No. of particles, Volume, and Temperature) ensemble for 1ns. (3) Subsequently, a 50,000,000-step (100 ns) production run was performed at constant pressure (1 bar) and 310K temperature after 1 ns NPT (Wang J. et al., 2020). During the MD simulations, the Particle Mesh Ewald (PME) (Wang et al., 2016) method was applied to account for the long-range electrostatic interactions, and the time step was set to 2 fs in all MD simulations. The salt bridge were calculated using VMD software (Humphrey et al., 1996). Further analyses (RMSD, RMSF, SASA, distances, angles and hydrogen bonds) were performed based on the resulting trajectories by GROMACS tools (Pronk et al., 2013). RMSD shows the deviation from the minimized crystal structure while RMSF shows the deviation from the mean structure over a dynamic ensemble. SASA measures the surface area of a molecule that is accessible to the solvent, generally water, and compares over the course of a simulation to detect solvent exposure events and changes to the protein surface. The SHAKE (Kr?Utler et al., 2015) algorithm was employed to calculate the bonds involving hydrogen atoms, for the H-bond interaction analysis, the two limiting factors are adopted as follows: (1) donor-acceptor distance is ≤ 3.5 Å, and (2) donor-hydrogen-acceptor angle is ≥ 150°.

MM-PBSA Calculation

Binding free energies of RBDs with hACE2 and mAbs were calculated using g_mmpbsa program. For each binding complex, 200 configurations were taken at an interval of 100 ps from the last 20 ns simulations. The conditions of MM-PBSA calculations are set, including solute dielectric constant = 2, solvent dielectric constant = 80, reference or vacuum dielectric constant = 1 and salt concentration = 0.15. The g_mmpbsa was used to calculate the polar solvation energies, non-polar solvation energies and calculates the free energy of the complex. The general expression of the term is

$\mathrm{\Delta }{G}_{bind}=G_{complex}-(G_{receptor}+G_{ligand})(1)$

Where G_complex is the total free energy complex, and G_receptor and G_ligand are total free energies of the isolated receptor and ligand in solvent, respectively.

$G_x=E_{MM}-TS+G_{solvation}(2)$

Where x is the receptor or ligand or complex. E_MM is the average molecular mechanics’ potential energy in a vacuum. TS refers to the entropic contribution to the free energy in a vacuum where T and S denote the temperature and entropy. G_solvation is the free energy of solvation.

$E=E_{bonded}+E_{nonbonded}=E_{bonded}+(E_{vdW}+E_{elec})(3)$

Where E_bonded is bonded interactions consisting of bond, angle, dihedral, and improper interactions. E_nonbonded includes both electrostatic and van der Waals interactions, which are depicted using a Coulomb and Lennard-Jones potential function, respectively. In addition, the free energy of solvation has been calculated including polar and nonpolar solvation energies, it that can be depicted as

$G_{solvation}=G_{polar}+G_{nonpolar}(3)$

Where G_polar and G_nonpolar are the electrostatic and non-electrostatic contributions to the solvation free energy, respectively. To check the statistical significance between the two complexes, P-value was calculated using t-test.

Results

The Diversity of Mutations in SARS-CoV-2 Whole-Genome Sequence

Our trajectory analysis of SARS-CoV-2 mutations in the COVID-19 pandemic was established on 64039 SARS-CoV-2 genome sequences (July 14, 2020) which were downloaded from the Global Initiative for Sharing All Influenza Data (GISAID) database. The SARS-CoV-2 sequences were aligned to the Wuhan reference genome with MAFFT (Katoh and Standley, 2013), then the amino acid changes were identified based on the sequence alignment. In general, most mutations are significantly concentrated in European and North American populations (Figure 1A). When observing the frequency of mutations in all genes, the amino acid changes are mostly distributed in ORF1a, ORF1b, N, and S proteins. Interestingly, the ORF3a protein has a superior mutation rate in North American and Oceanian populations accounting for 14.81 and 10.22% of the total protein, respectively (Figure 1A). The ORF3a protein is also called “accessory protein,” it can convert the environment inside the infected cell and make holes on the infected cell membrane.

FIGURE 1

Figure 1. The distribution of Missense Mutations in SARS-CoV-2. A stacked histogram shows all missense mutations of SARS-CoV-2 from different continents of the world. (B) A chord diagram shows the missense mutation shared between different proteins of SARS-CoV-2, and the color scheme is the same as that in (A). (C) A topology and cartoon representation of SARS-CoV-2 homo-trimeric spike (S) glycoprotein (PDB ID: 6VSB). NTD, N-terminal domain; FP, fusion peptide; HR1, heptad repeat 1; CH, central helix; CD, connector domain; HR2, heptad repeat 1; TM, transmembrane region; CT, cytoplasmic tail; S1/S2, protease cleavage sites; Some important mutation sites and residue intervals were marked on the figure. (D) A Venn diagram of main amino acid mutation sites of the S protein, the number of mutations, and co-mutations were shown.

The co-mutation may affect the cooperation between the various proteins, which may lead to meaningful virus evolution and pose challenges to the development of antibodies. To investigate the co-mutation in SARS-CoV-2, we calculated the co-mutation in all proteins and found most mutants were not single-site mutations (Figure 1B). Since SARS-CoV-2 infects humans through the binding of trimeric spike glycoprotein (PDB ID: 6VSB) to hACE2, then we counted all amino acid variants on S protein (Figure 1C and Supplementary Table 1). The variant D614G accounts for 75.92% of 64039 SARS-CoV-2 (Korber et al., 2020), following by D936Y accounts for 1.11%. And more notably, N439K is the highest frequency of variant in the RBD region, which accounts for 0.72% of all SARS-CoV-2. Furthermore, there are only 16 mutations on S protein with a mutation rate of more than 0.16% (Supplementary Table 1), and those mutants occur frequently together with the variant D614G (Figure 1D). All the N439K are included in the D614G interestingly, which is mainly concentrated in Europe since March 2020. To conclude, the amino acid changes seem to have been accumulated progressively over time, among which N439K is the most dominant variant in the RBD region and should be preferentially employed to characterize the influence of mutations on pathogen evolution.

The N439K-Mutated RBD Binds hACE2 With Higher Affinity Than Wild Type

Crystal structures show that the high contact area of the SARS-CoV-2 RBD-hACE2 complex is mainly concentrated in the RBM region (Shang et al., 2020; Wang Q. et al., 2020), and denoted as CR1, CR2, and CR3 (Figure 2A; Wang Y. et al., 2020). The PROCHECK analysis show both wild-type and N439K S protein model are in good quality, the detail information can be found in Supplementary File 2. To further verify the convergence of MD simulations equilibrium, we estimated the root mean square deviations (RMSD) of backbone atoms relative to the corresponding crystal structure, and the wild complex had a relatively smaller average RMSD than the N439K-mutated complex (2.2 vs. 2.4 Å) (Figure 2B and Table 1). The structural compactness of each system was elucidated by estimating the radius of gyration (R_g), the average R_gvalues were similar for all the same systems (Table 1). We also calculated the solvent accessible surface area (SASA) which indicates the solvent exposure degree, and SASA values of wild and mutant complexes were 35716 Ų and 36121 Ų, respectively (Table 1 and Supplementary Table 3). Although the flexibility patterns of residues in the N439K-mutated RBD-hACE2 complex display similar fluctuations with wild-type complex (Figure 2C), certain regions of the two complexes show differences in flexibility.

FIGURE 2

Figure 2. The RBD-hACE2 interaction profile of MD simulations. (A) The structure of SARS-CoV-2 RBD-hACE2 (PDB ID: 6M0J). SARS-CoV-2 RBD core is shown in deepsalmon and its interacting hACE2 is colored greencyan, and RBM of RBD is represented in yellow. Receptor-ligand interface includes the N-terminal (CR1), the central region (CR2) and the binding loop (CR3), ASN439 (yellow) to LYS439 (red) in CR3. (B) The RMSDs of the backbone atoms of both RBD-hACE2 complexes, the RBD-hACE2 is colored orange and RBD(N439K)-hACE2 is colored blue. (C) The RMSFs of C_αatoms of both RBD-hACE2 complexes, where RBD-hACE2 is colored orange and RBD(N439K)-hACE2 is colored blue. (D,E) Hydrogen bonds and salt bridges are shown as dotted lines, RBD (deepsalmon), and hACE2 (greencyan) residues are described as sticks, respectively. The orange box indicates wild-type RBD-hACE2 and mutant-type RBD-hACE2 is a blue box.

TABLE 1

Table 1. The average RMSD, solvent accessible surface area (SASA), and radius of gyration for the simulated systems.

Subsequently, we performed the hydrogen bond (H-bond) and salt bridge analyses based on the trajectories in MD simulations. It was observed that the N439K-mutated RBD-hACE2 complex (9.68 ± 3.24) can form more hydrogen bonds than the wild-type (7.26 ± 2.59) during the MD trajectory (P = 2.43 × 10^–138) (Figure 2D and Supplementary Table 4). Hydrogen bonds from the N-terminal (CR1) to the central region (CR2) of the interface were almost the same during MD simulations, and three additional main-chain hydrogen bonds forms at Gly496, Gln498, and Val503 in the binding loop (CR3), causing the ridge to take a more compact conformation and the loop to move closer to hACE2 (Figure 2D). At the SARS-CoV-2 RBD-hACE2 interface, a strong salt bridge (Pylaeva et al., 2018) between Lys417 of the RBD and Asp30 of hACE2 had been confirmed (Lan et al., 2020; Figure 2E). The mutation of N439K in the RBD formed a new salt bridge with Glu329 of hACE2 (3.6 Å), and a weak salt bridge between Lys439 and Asp442 of SARS-CoV-2 RBD (Figure 2D). Burial of these salt bridges in hydrophobic environments on virus binding would enhance their energy owing to a reduction in the dielectric constant.

To estimate the stabilization of binding systems, We also calculated the molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) in both wild and mutant SARS-CoV-2 RBD-hACE2 complexes (Kumari et al., 2014). Overall, the binding free energy (ΔG_bind) decomposition analysis divulged into various free energies (Figure 3B and Supplementary Table 5). The binding energy ΔG_bind (−1526.17 ± 133.13 kj/mol) of the N439K-hACE2 was higher in magnitude as compared to the wild-type RBD-hACE2 (−1084.06 ± 80.23kj/mol) (Figure 3A), in which the electrostatic energy between wild and mutant types had a significant difference, which the average ΔE_elec (−1876.25 ± 45.23 kj/mol) of the N439K-hACE2 also higher than that of wild-type (−1301 ± 23.14 kj/mol) (Figure 2B). To test the reliability of our simulation, we increase the MD simulation time to 200 ns, the binding energy of 200 ns MD simulation show a similar result, the ΔG_bind (−1597.25 ± 179.57 kj/mol) of the N439K-hACE2 also higher than that of wild-type (−959.29 ± 130.36 kj/mol) (Supplementary Table 5). Meanwhile, we calculated binding energy between RBM (residue 438–502) and hACE2 using MM-PBSA, and binding energies of the two complexes were (−860.94 ± 99.45 kj/mol and (−398.06 ± 111.76 kj/mol, respectively (Figure 3C and Supplementary Table 6). Comparing the binding free energy of RBD-hACE2 and RBM-hACE2, it should be noted that the changes in energy are mainly concentrated in the RBM-hACE2 region.

FIGURE 3

Figure 3. Energy components for binding free energy of both wild and mutated RBD/RBM-hACE2 complexes. (A) The binding free energies for SARS-CoV-2 RBD-hACE2 (including wild type and variantN439K) using MM-PBSA. Total binding energy (ΔGbind). (B) Energy components for the binding free energy of RBD-hACE2. The intermolecular van der Waals (ΔEvdw); Electrostatic interactions (ΔEelec); Polar solvation energy (ΔGpol) and apolar (non-polar) solvation energy (ΔGapol). (C) The binding free energies for SARS-CoV-2 RBM-hACE2 (including wild type and variant N439K), using orange for wild-type and blue for mutant type. (D) Decomposition of ΔGbind into contributions from individual residues (438–502) of RBM before and after mutation. All units are reported in kj/mol.

Subsequently, we explored the critical residues involved in the RBM-hACE2 binding by performing the per-residue decomposition of binding free energy (Figure 3D), and the contribution energy of per-residue to the total binding energy was compared before and after mutation. It is evident from the line graph (Figure 3D) that Lys439 hotspot residue show more contribution to the binding free energy than ASN439 and its contribution energy changed from −1.93 to −249.85 kj/mol (Figure 3D). The total binding free energy change between wild-type and RBD-hACE2 (N439K) complexes is 462.88 kj/mol (Figure 3C), mostly concentrated in the electrostatic energy (Figure 3B). Taken together, comparing the interaction interfaces of the N439K-mutated and wild RBD-hACE2 complexes reveals the change from ASN439 to LYS439 might result in a tighter association because of the new salt bridge formation and higher affinity.

N439K Became Resistant to SARS-CoV-2 Neutralizing Antibody REGN10987

To investigate the antigenicity of the N439K mutant, we exploited 100 ns MD simulations of the binary complexes of hACE2 with neutralizing monoclonal antibodies REGN10987 and CB6 complexes. Specific neutralizing mAbs can prevent the virus from binding to hACE2 by neutralizing SARS-CoV-2. Although most of these highly binding sites are in RBM, the two neutralizing antibodies show different high contact regions with RBD (Figure 4A).

FIGURE 4

Figure 4. Structural and energetic details of both wild and mutant RBD-mAbs interactions. (A) Crystal structures of RDB-CB6/REGN10987 complexes, the RBD is colored red, CB6 heavy and light chains are represented as marine and green, respectively, REGN10987 heavy chain is colored yellow and light is blue, and the 439 residues are described as the sphere. (B) Characteristic dynamic fluctuations of both RBD-REGN10987 and RBD(N439K)-REGN10987 complexes. Mutant-type (100 ns) and wild-type (100 ns) are colored by orange and blue, respectively. (C) Dynamic conformations are projected on to the principal vectors (PC1 and PC2). Red and blue indicate mutant-type and wild-type 100 ns MD trajectories, respectively. (D) The RMSD of the receptor-binding motif in four complexes during the 100-ns MD simulations. (E) The binding free energies for both complexes of the mAb REGN10987 (including heavy and light chains), the color schemes are the same in (A,C). (F) The binding free energies of 200 configurations at an interval of 100 ps from the last 20 ns simulations. The t-test was conducted to check the statistical significance of the difference between two systems of binding free energies. A p-value of <0.05 indicates that the difference is statistically significant (95% confidence interval). The color scheme is the same as that in (C).

After the 100 ns MD simulation, the final conformational distribution indicated the configuration fluctuations of RBD were mainly concentrated in the CR3, whereas the REGN10987 antibody was relatively stable (Figure 4B). Overall, SARS-CoV-2 RBD and REGN10987 undergo symmetric twist and antisymmetric hinge-bending motions about the axis of the N-terminal helix, corresponding to the lowest quasi-harmonic models, PC1, and PC2, from principal component analysis (PCA) (Ichiye and Karplus, 1991). Although the overall conformational dynamics were almost the same between the wild and mutant complexes, the average conformation of the variant structure has changed relative to wild complex when the dynamic configurations are projected onto the two principal vectors (Figure 4C). From results of RMSD, the wild and mutant systems had reached a stable state from their respective MD trajectories, CB6 systems had lower average RMSD values than the REGN10987 complexes, and a higher RMSD was calculated in N439K-REGN10987 than wild-type (Figure 4D). The MD simulation results of the CB6-RBD complex can be found in Supplementary Table 7.

To further explain the ability of heavy and light chains of mAbs to neutralize viruses, we had taken 200 structures from the stable region of the last 20ns trajectories to calculate binding free energy with RBD, respectively. It revealed that the estimated binding free energy of wild-type CB6-RBD complexes (−129.73 kj/mol) is higher than N439K-mutated CB6-RBD complexes (−51.28 kj/mol) in the heavy chain, the unfavorable contribution from ΔG_pol (491.56 kj/mol) was relatively lower compared to mutant type (703.82 kj/mol) (Table 2). Overall, our research suggested that the N439K reduced the sensitivity to the CB6 mAb. Subsequently, we have investigated the binding free energy of REGN10987-RBD complexes, the ΔG_bind of N439K-mutated REGN10987-RBD complexes (24.10 kj/mol) was found to be lower (P = 2.22e-16) than wild-type (−85.65 kj/mol) in the heavy chain (Figure 4F and Table 2). Similarly, the light chain had a little difference in binding free energy with RBD, but both systems could not form an effective combination (Figure 4E and Table 2). Furthermore, the present 100 ns trajectory revealed that polar solvation and electrostatic interaction might be the main factors to lose the ability of the REGN10987 antibody to neutralize COVID-19 (Table 2).

TABLE 2

Table 2. Energetic components of binding energy for SARS-CoV-2 RBD-mAbs complexes (kj/mol).

Discussion

The outbreak of COVID-19 has caused a severe strain on the public health system in many countries. The SARS-CoV-2 virus is expected to continue evolving in human populations. Close monitoring of circulating virus strains is of unquestionable importance to inform research and development of antibodies and therapeutics. Herein, we have studied the mechanism of binding of RBD-hACE2 and RBD-mAbs by using an atomistic molecular dynamics simulation in conjunction with molecular mechanics/Poisson-Boltzmann surface area scheme. Our study shows that the N439K influences the affinity of both RBD-hACE2 and RBD-mAbs complexes which is favored by the intermolecular van der Waals, electrostatic interactions, and polar solvation free energy. Our findings facilitate the development of decoy ligands and neutralizing antibodies for suppression of viral infection.

MD simulations and MM-PBSA results showed that the binding ability of the N439K-mutated RBD with hACE2 was enhanced. In other words, N439K increases spike affinity for ACE2 significantly, which is consistent with the surface plasmon resonance (SPR) result (Thomson et al., 2021). The replacement of ARG426 (SARS-CoV RBD) with ASN439 (SARS-CoV-2 RBD) appeared to weaken the interaction by eliminating one important salt bridge with ASP329 on hACE2 and reduced the affinity (Lan et al., 2020; Yan et al., 2020). The N439K (ASN439 to LYS439) SARS-CoV-2 variant form a new salt bridge at the RBD-hACE2 interface through simulation calculation, which is consistent with the previous study (Thomson et al., 2021). Thomson used SPR to evaluate binding of recombinant N439K RBD protein to recombinant hACE2, indicating that acquisition of the N439K mutation enhances hACE2 binding. Meanwhile, the reduced binding of REGN10987 mAb to the variant N439K RBD, which was also confirmed by bio-layer interferometry analysis (Thomson et al., 2021). Compared with traditional experimental methods, MD simulations require less time, and can explain structure and energy changes at molecular level, such as salt bridge and binging free energy. Further, the highly consistent results with experimental methods verify the power of our method. As the SARS-CoV-2 virus is expected to continue evolving in populations, molecular dynamics simulation is characteristic of high speed and low cost, which is especially suitable for high throughput screening of high-risk mutation in S protein.

A few other circulating RBD mutations have become prominent since N439K first emerged. Independent lineages of SARS-CoV-2 have recently been reported: UK-B.1.1.7, South Africa-B.1.351 and Brail-P.1. The B.1.1.7 variants with increased transmission have 9 amino-acid changes in Spike, including N501Y, and N501Y compromises neutralization by many antibodies with public V-region IGHV3-53 (Supasa et al., 2021). The B.1.351 variants of SARS-CoV-2 containing multiple mutations in Spike are now dominant in South Africa, indicating the potential for impaired efficacy of potential monoclonal antibodies and vaccines (Li et al., 2021). Among the E484K, K417N, and N501Y mutations in the receptor-binding domain of Spike caused widespread escape from monoclonal antibodies (Zhou et al., 2021). The Y435F mutation provides evidence of animal-to-human transmission in mink farms (Oude Munnink et al., 2021). However, these variants with high rates of infection have gained notice, further are studied the impact of own monoclonal antibodies and vaccines. Due to widely spread, it’s difficult to play an early warning role. Whether it is possible to find a way to quickly and effectively monitor the spread of SARS-CoV-2 is still an urgent problem.

In general, vaccinations are being deployed worldwide, which can widely reduce the spread of SARS-CoV-2 among the population. However, as the duration of a SARS-CoV-2 pandemic has prolonged, there have been reported multiple variants around the world. Here, we describe a circulating RBD mutation N439K that maintains a high affinity with hACE2 while evading antibody-mediated immune response using bioinformatics method. More similar reports for SARS-CoV-2 variants have been pointed, and as more and more people develop immune responses against virus via infection or vaccination, the monitoring of SARS-CoV-2 escape mutants is crucial, and may require new vaccine preparations that address the variants circulating globally.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found here: Publicly available datasets were analyzed in this study. The original sequencing of SARS-CoV-2 can be downloaded from downloaded from the GISAID (https://www.gisaid.org/) database, and the wild-type crystal structure of spike protein binding to ACE2 and mAbs were downloaded from the Protein Data Bank (PDB ID: 6M0J, 6XDG, and 7C01).

Author Contributions

QJ and HN conceived the project. WZ and CX collected genome sequences and performed multiple sequence alignment. WZ, QJ, CX, PW, ZX, RC, XJ, GX, YG, GX, AA, and LJ contributed to molecular dynamics and binding free energy simulations. QJ, HN, WZ, and CX wrote the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was funded by the National Natural Science Foundation of China (Nos. 61822108, 62032007 and 62041102 to QJ) and the Emergency Research Project for COVID-19 of Harbin Institute of Technology (No. 2020-001 to QJ).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.697035/full#supplementary-material

Supplementary Table 1 | Statistics of the mutations of SARS-CoV-2 from GISAID database.

Supplementary Table 3 | Details on both RBD-hACE2 complexes for all simulated systems.

Supplementary Table 4 | The number of hydrogen bonds between RBD and hACE2 protein in N439K and wild-type RBD-hACE2 complex.

Supplementary Table 5 | Summary on energetic components of binding energy for SARS-CoV-2 RBD-ACE2 complexes.

Supplementary Table 6 | Details on energetic components of binding free energy for RBD/RBM-hACE2 complexes.

Supplementary Table 7 | Details on binding free energy components for RBD-mAbs complexes.

Supplementary File 2 | The quality control results of both wild-type and N439K S protein model.

Footnotes

1. ^ https://covid19.who.int/
2. ^ https://www.rcsb.org/

References

Asmal, M., Hellmann, I., Liu, W., Keele, B. F., Perelson, A. S., Bhattacharya, T., et al. (2011). A signature in HIV-1 envelope leader peptide associated with transition from acute to chronic infection impacts envelope processing and infectivity. PloS One 6:e23673. doi: 10.1371/journal.pone.0023673

PubMed Abstract | CrossRef Full Text | Google Scholar

Bowie, J., Luthy, R., and Eisenberg, D. (1991). A method to identify protein sequences that fold into a known three-dimensional structure. Science 253, 164–170.

Google Scholar

Bricault, C. A., Yusim, K., Seaman, M. S., Yoon, H., Theiler, J., Giorgi, E. E., et al. (2019). HIV-1 neutralizing antibody signatures and application to epitope-targeted vaccine design. Cell Host Microbe 25, 5972.e8. doi: 10.1016/j.chom.2018.12.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, Y., Su, B., Guo, X., Sun, W., Deng, Y., Bao, L., et al. (2020). Potent neutralizing antibodies against SARS-CoV-2 identified by high-throughput single-cell sequencing of convalescent patients’ B cells. Cell 182, 73–84.e16. doi: 10.1016/j.cell.2020.05.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Hansen, J., Baum, A., Pascal, K. E., Russo, V., Giordano, S., Wloga, E., et al. (2020). Studies in humanized mice and convalescent humans yield a SARS-CoV-2 antibody cocktail. Science 369, 1010–1014. doi: 10.1126/science.abd0827

PubMed Abstract | CrossRef Full Text | Google Scholar

Hoffmann, M., Kleine-Weber, H., Schroeder, S., Krüger, N., Herrler, T., Erichsen, S., et al. (2020). SARS-CoV-2 cell entry depends on ACE2 and TMPRSS2 and is blocked by a clinically proven protease inhibitor. Cell 181, 271–280.e8. doi: 10.1016/j.cell.2020.02.052

PubMed Abstract | CrossRef Full Text | Google Scholar

Humphrey, W., Dalke, A., and Schulten, K. (1996). VMD: visual molecular dynamics. J. Mol. Graph. 14, 33–38.

Google Scholar

Ichiye, T., and Karplus, M. (1991). Collective motions in proteins: a covariance analysis of atomic fluctuations in molecular dynamics and normal mode simulations. Proteins 11, 205–217.

Google Scholar

Ju, B., Zhang, Q., Ge, J., Wang, R., Sun, J., Ge, X., et al. (2020). Human neutralizing antibodies elicited by SARS-CoV-2 infection. Nature 584, 115–119. doi: 10.1038/s41586-020-2380-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Katoh, K., and Standley, D. M. (2013). MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol. Biol. Evol. 30, 772–780. doi: 10.1093/molbev/mst010

PubMed Abstract | CrossRef Full Text | Google Scholar

Korber, B., Fischer, W. M., Gnanakaran, S., Yoon, H., Theiler, J., Abfalterer, W., et al. (2020). Tracking changes in SARS-CoV-2 spike: evidence that D614G increases infectivity of the COVID-19 Virus. Cell 182, 812–827.e19. doi: 10.1016/j.cell.2020.06.043

PubMed Abstract | CrossRef Full Text | Google Scholar

Kr?Utler, V., Gunsteren, W. F. V., and Hünenberger, P. H. (2015). A fast SHAKE algorithm to solve distance constraint equations for small molecules in molecular dynamics simulations. J. Comput. Chem. 22, 501–508.

Google Scholar

Kumari, R., Kumar, R., and Lynn, A. (2014). g_mmpbsa–a GROMACS tool for high-throughput MM-PBSA calculations. J. Chem. Inf. Model. 54, 1951–1962.

Google Scholar

LaBranche, C. C., Henderson, R., Hsu, A., Behrens, S., Chen, X., Zhou, T., et al. (2019). Neutralization-guided design of HIV-1 envelope trimers with high affinity for the unmutated common ancestor of CH235 lineage CD4bs broadly neutralizing antibodies. PLoS Pathog. 15:e1008026. doi: 10.1371/journal.ppat.1008026

PubMed Abstract | CrossRef Full Text | Google Scholar

Lan, J., Ge, J., Yu, J., Shan, S., Zhou, H., Fan, S., et al. (2020). Structure of the SARS-CoV-2 spike receptor-binding domain bound to the ACE2 receptor. Nature 581, 215–220. doi: 10.1038/s41586-020-2180-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Laskowski, A. R., Macarthur, W. M., Moss, S. D., and Thornton, M. J. (1993). PROCHECK: a program to check the stereochemical quality of protein structures. J. Appl. Cryst. 26, 283–291.

Google Scholar

Li, F., Luo, M., Zhou, W., Li, J., Jin, X., Xu, Z., et al. (2020). Single cell RNA and immune repertoire profiling of COVID-19 patients reveal novel neutralizing antibody. Protein Cell [Epub ahead of print], 1–5. doi: 10.1007/s13238-020-00807-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Q., Nie, J., Wu, J., Zhang, L., Ding, R., Wang, H., et al. (2021). No higher infectivity but immune escape of SARS-CoV-2 501Y.V2 variants. Cell 184, 2362–2371.e9. doi: 10.1016/j.cell.2021.02.042

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Q., Wu, J., Nie, J., Zhang, L., Hao, H., Liu, S., et al. (2020). The Impact of Mutations in SARS-CoV-2 Spike on Viral Infectivity and Antigenicity. Cell 182, 1284.e–1294.e. doi: 10.1016/j.cell.2020.07.012 1284-1294.e9

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, W., Moore, M. J., Vasilieva, N., Sui, J., Wong, S. K., Berne, M. A., et al. (2003). Angiotensin-converting enzyme 2 is a functional receptor for the SARS coronavirus. Nature 426, 450–454.

Google Scholar

Lilkova, E., Petkov, P., Ilieva, N., and Litov, L. (2015). The PyMOL Molecular Graphics System Version 2.0.

Google Scholar

Liu, L., Wang, P., Nair, M. S., Yu, J., Rapp, M., Wang, Q., et al. (2020). otent neutralizing antibodies against multiple epitopes on SARS-CoV-2 spike. Nature 584, 450–456. doi: 10.1038/s41586-020-2571-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, R., Zhao, X., Li, J., Niu, P., Yang, B., Wu, H., et al. (2020). Genomic characterisation and epidemiology of 2019 novel coronavirus: implications for virus origins and receptor binding. Lancet 395, 565–574. doi: 10.1016/S0140-6736(20)30251-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Melo, F., Devos, D., Depiereux, E., and Feytmans, E. (1997). ANOLEA: a www server to assess protein structures. Proc. Int. Conf. Intell. Syst. Mol. Biol. 5, 187–190.

Google Scholar

Morse, J. S., Lalonde, T., Xu, S., and Liu, W. R. (2020). Learning from the past: possible urgent prevention and treatment options for severe acute respiratory infections caused by 2019-nCoV. Chembiochem 21, 730–738. doi: 10.1002/cbic.202000047

PubMed Abstract | CrossRef Full Text | Google Scholar

Oude Munnink, B. B., Sikkema, R. S., Nieuwenhuijse, D. F., Molenaar, R. J., Munger, E., Molenkamp, R., et al. (2021). Transmission of SARS-CoV-2 on mink farms between humans and mink and back to humans. Science 371, 172–177. doi: 10.1126/science.abe5901

PubMed Abstract | CrossRef Full Text | Google Scholar

Pontius, J., Richelle, J., and Wodak, S. J. (1996). Deviations from standard atomic volumes as a quality measure for protein crystal structures. J. Mol. Biol. 264, 121–136.

Google Scholar

Price, D. J., and Brooks, C. L. (2004). A modified TIP3P water potential for simulation with Ewald summation. J. Chem. Phys. 121, 10096–10103.

Google Scholar

Pronk, S., Páll, S., Schulz, R., Larsson, P., Bjelkmar, P., Apostolov, R., et al. (2013). GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics 29, 845–854.

Google Scholar

Pylaeva, S., Brehm, M., and Sebastiani, D. (2018). Salt bridge in aqueous solution: strong structural motifs but weak enthalpic effect. Sci. Rep. 8:13626. doi: 10.1038/s41598-018-31935-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Saha, P., Banerjee, A. K., Tripathi, P. P., Srivastava, A. K., and Ray, U. (2020). A virus that has gone viral: amino acid mutation in S protein of Indian isolate of coronavirus COVID-19 might impact receptor binding, and thus, infectivity. Biosci. Rep. 40:BSR20201312. doi: 10.1042/BSR20201312

PubMed Abstract | CrossRef Full Text | Google Scholar

Shang, J., Ye, G., Shi, K., Wan, Y., Luo, C., Aihara, H., et al. (2020). Structural basis of receptor recognition by SARS-CoV-2. Nature 581, 221–224. doi: 10.1038/s41586-020-2179-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Sheikh, J. A., Singh, J., Singh, H., Jamal, S., Khubaib, M., Kohli, S., et al. (2020). Emerging genetic diversity among clinical isolates of SARS-CoV-2: lessons for today. Infect. Genet. Evol. 84:104330. doi: 10.1016/j.meegid.2020.104330

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, R., Shan, C., Duan, X., Chen, Z., Liu, P., Song, J., et al. (2020). A human neutralizing antibody targets the receptor-binding site of SARS-CoV-2. Nature 584, 120–124. doi: 10.1038/s41586-020-2381-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Sui, J., Aird, D. R., Tamin, A., Murakami, A., Yan, M., Yammanuru, A., et al. (2008). Broadening of neutralization activity to directly block a dominant antibody-driven SARS-coronavirus evolution pathway. PLoS Pathog. 4:e1000197. doi: 10.1371/journal.ppat.1000197

PubMed Abstract | CrossRef Full Text | Google Scholar

Supasa, P., Zhou, D., Dejnirattisai, W., Liu, C., Mentzer, A. J., Ginn, H. M., et al. (2021). Reduced neutralization of SARS-CoV-2 B.1.1.7 variant by convalescent and vaccine sera. Cell 184, 2201–2211.e7. doi: 10.1016/j.cell.2021.02.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Tang, X.-C., Agnihothram, S. S., Jiao, Y., Stanhope, J., Graham, R. L., Peterson, E. C., et al. (2014). Identification of human neutralizing antibodies against MERS-CoV and their role in virus adaptive evolution. Proc. Natl. Acad. Sci. U. S. A. 111, E2018–E2026. doi: 10.1073/pnas.1402074111

PubMed Abstract | CrossRef Full Text | Google Scholar

ter Meulen, J., van den Brink, E. N., Poon, L. L. M., Marissen, W. E., Leung, C. S. W., Cox, F., et al. (2006). Human monoclonal antibody combination against SARS coronavirus: synergy and coverage of escape mutants. PLoS Med. 3:e237. doi: 10.1371/journal.pmed.0030237

PubMed Abstract | CrossRef Full Text | Google Scholar

Thomson, E. C., Rosen, L. E., Shepherd, J. G., Spreafico, R., Filipe, A. D. S., Wojcechowskyj, J. A., et al. (2021). Circulating SARS-CoV-2 spike N439K variants maintain fitness while evading antibody-mediated immunity. Cell 184, 1171–1187.e20.

Google Scholar

Tortorici, M. A., and Veesler, D. (2019). Structural insights into coronavirus entry. Adv. Virus Res. 105, 93–116. doi: 10.1016/bs.aivir.2019.08.002

PubMed Abstract | CrossRef Full Text | Google Scholar

van Dorp, L., Acman, M., Richard, D., Shaw, L. P., Ford, C. E., Ormond, L., et al. (2020). Emergence of genomic diversity and recurrent mutations in SARS-CoV-2. Infect. Genet. Evol. 83:104351. doi: 10.1016/j.meegid.2020.104351

PubMed Abstract | CrossRef Full Text | Google Scholar

Walls, A. C., Park, Y.-J., Tortorici, M. A., Wall, A., McGuire, A. T., and Veesler, D. (2020). Structure, function, and antigenicity of the SARS-CoV-2 spike glycoprotein. Cell 181, 281–292.e6. doi: 10.1016/j.cell.2020.02.058

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H., Gao, X., and Fang, J. (2016). Multiple staggered mesh ewald: boosting the accuracy of the smooth particle mesh Ewald method. J. Chem. Theor. Comput. 12, 5596–5608.

Google Scholar

Wang, J., Xu, X., Zhou, X., Chen, P., and Hao, P. (2020). Molecular simulation of SARS-CoV-2 spike protein binding to pangolin ACE2 or human ACE2 natural variants reveals altered susceptibility to infection. J. Gen. Virol. 101, 921–924.

Google Scholar

Wang, P., Jin, X., Zhou, W., Luo, M., Xu, Z., Xu, C., et al. (2021). Comprehensive analysis of TCR repertoire in COVID-19 using single cell sequencing. Genomics 113, 456–462. doi: 10.1016/j.ygeno.2020.12.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Q., Zhang, Y., Wu, L., Niu, S., Song, C., Zhang, Z., et al. (2020). Structural and functional basis of SARS-CoV-2 entry by using human ACE2. Cell 181, 894–904.e9. doi: 10.1016/j.cell.2020.03.045

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Liu, M., and Gao, J. (2020). Enhanced receptor binding of SARS-CoV-2 through networks of hydrogen-bonding and hydrophobic interactions. Proc. Natl. Acad. Sci. U. S. A. 117, 13967–13974. doi: 10.1073/pnas.2008209117

PubMed Abstract | CrossRef Full Text | Google Scholar

Waterhouse, A., Bertoni, M., Bienert, S., Studer, G., Tauriello, G., Gumienny, R., et al. (2018). SWISS-MODEL: homology modelling of protein structures and complexes. Nucleic Acids Res. 46, W296–W303. doi: 10.1093/nar/gky427

PubMed Abstract | CrossRef Full Text | Google Scholar

Wrapp, D., Wang, N., Corbett, K. S., Goldsmith, J. A., Hsieh, C.-L., Abiona, O., et al. (2020). Cryo-EM structure of the 2019-nCoV spike in the prefusion conformation. Science 367, 1260–1263. doi: 10.1126/science.abb2507

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, F., Zhao, S., Yu, B., Chen, Y. M., Wang, W., Song, Z. G., et al. (2020). A new coronavirus associated with human respiratory disease in China. Nature 579, 265–269. doi: 10.1038/s41586-020-2008-3

CrossRef Full Text | Google Scholar

Yan, R., Zhang, Y., Li, Y., Xia, L., Guo, Y., and Zhou, Q. (2020). Structural basis for the recognition of SARS-CoV-2 by full-length human ACE2. Science 367, 1444–1448. doi: 10.1126/science.abb2762

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, D., Dejnirattisai, W., Supasa, P., Liu, C., Mentzer, A. J., Ginn, H. M., et al. (2021). Evidence of escape of SARS-CoV-2 variant B.1.351 from natural and vaccine induced sera. Cell 184, 2348–2361.e6. doi: 10.1016/j.cell.2021.02.037

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, P., Wang, H., Fang, M., Li, Y., Wang, H., Shi, S., et al. (2019). Broadly resistant HIV-1 against CD4-binding site neutralizing antibodies. PLoS Pathog. 15:e1007819. doi: 10.1371/journal.ppat.1007819

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, P., Yang, X.-L., Wang, X.-G., Hu, B., Zhang, L., Zhang, W., et al. (2020). A pneumonia outbreak associated with a new coronavirus of probable bat origin. Nature 579, 270–273. doi: 10.1038/s41586-020-2012-7

PubMed Abstract | CrossRef Full Text | Google Scholar