Allosteric Cross-Talk among Spike's Receptor-Binding Domain Mutations of the SARS-CoV-2 South African Variant Triggers an Effective Hijacking of Human Cell Receptor.

The rapid and relentless emergence of novel highly transmissible SARS-CoV-2 variants, possibly decreasing vaccine efficacy, currently represents a formidable medical and societal challenge. These variants frequently hold mutations on the Spike protein's receptor-binding domain (RBD), which, binding to the angiotensin-converting enzyme 2 (ACE2) receptor, mediates viral entry into host cells. Here, all-atom molecular dynamics simulations and dynamical network theory of the wild-type and mutant RBD/ACE2 adducts disclose that while the N501Y mutation (UK variant) enhances the Spike's binding affinity toward ACE2, the concomitant N501Y, E484K, and K417N mutations (South African variant) aptly adapt to increase SARS-CoV-2 propagation via a two-pronged strategy: (i) effectively grasping ACE2 through an allosteric signaling between pivotal RBD structural elements and (ii) impairing the binding of antibodies elicited by infected or vaccinated patients. This information unlocks the molecular terms and evolutionary strategies underlying the increased virulence of emerging SARS-CoV-2 variants, setting the basis for developing the next-generation anti-COVID-19 therapeutics.

The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), the etiological agent of coronavirus disease 19 (COVID-19), has infected as of April 30, 2021 about 150 million patients, causing over 3 million deaths worldwide. Owing to an unprecedentedly intense and relentless scientific effort, a variety of vaccines and monoclonal antibodies are becoming available for COVID-19 prophylaxis and therapeutic treatment.[1−3]

Similar to other β-coronaviruses (β-CoVs), the receptor-binding domain (RBD) of the homotrimeric viral spike (S) protein of SARS-CoV-2 mediates the molecular recognition and the binding to the human cellular receptor, angiotensin-converting enzyme 2 (ACE2),[4,5] thus triggering SARS-CoV-2 entry into host cells. As such, the S-protein has been the object of burgeoning research interest, becoming the prominent target for antibody development. This prompted an exhaustive experimental[6−8] and computational[9−18] assessment of the molecular interactions between the S-protein and ACE2.

The worldwide continuous and uncontrolled transmission of SARS-CoV-2 set the condition for its rapid evolution into more infectious variants. As an example, one of the first S-protein mutations, D614G, characterized by an enhanced transmissibility, has rapidly become dominant.[19] As well, other alarming strains have emerged in United Kingdom (lineage B.1.1.7),[20] South Africa (lineage B.1.351),[21] and Brazil (lineage P.1),[22] hereafter termed the UK, SA, and BR variants, respectively. Ultimately, a new dire Indian variant (lineage B.1.617) came to the fore. These lineages are the object of rising concerns owing to their increased transmissibility and/or their potential ability to elude infection- or vaccine-induced immunity.

As concerns the most prominent nonsynonymous mutations placed in the S-protein’s RBD, most SARS-CoV-2 variants share the N501Y substitution (Figure ), most likely implicated into an enhanced binding affinity toward ACE2,[23−25] although preliminary reports indicate that this variant retains vaccine efficacy.[26] In addition to N501Y, the SA variant also exhibits the E484K and K417N RBD mutations. E484, the most frequently mutated residue in COVID-19 patients, becomes E484K in the SA and BR and E484Q in the Indian strains. As well, mutation of K417, either to N or a T, is shared by the SA and BR variants, respectively. These mutations have been linked to viral escape from mAbs developed by vaccinated or infected patients.[2,27,28]

Figure 1

Representative structures of the complex between the South African (SA) SARS-CoV-2 variant of the receptor-binding domain (RBD, pink, with the receptor-binding motif (RBM) highlighted in green) and the angiotensin-converting enzyme 2 (ACE2, blue) as extracted from molecular dynamics trajectories. The three N501Y, E484K, and K417N mutations sites are circled in red, yellow, and black, respectively. The insets show a comparison of the key intermolecular interactions at the mutation sites in the wild-type (WT) and SA RBD/ACE2 complexes with residues depicted in licorice and hydrogen bonds displayed as dashed lines.

Aiming to dissect at the atomic level the role of RBD mutations on the recognition of ACE2, we performed cumulative 15 μs all-atom molecular dynamics (MD) simulations of S-protein RBD/ACE2 complexes considering the RBD’s mutations present in the SA variant either concurrently or singularly.

Namely, we first built the adduct between ACE2 and RBD carrying the N501Y, E484K, and K417N substitutions of the SA lineage (hereafter termed SARBD/ACE2). Next, to inspect the role of each mutation, we built three distinct RBD/ACE2 models carrying N501Y (N501YRBD/ACE2 or UKRBD/ACE2), E484K (E484KRBD/ACE2), and K417N (K417NRBD/ACE2), ultimately comparing them with the WT RBD/ACE2 adduct (hereafter named RBD/ACE2). As a result, all the systems retain stable interactions at the RBD/ACE2 interface when performing 2.5 μs MD simulations for each system (Figure S1). Most of the RBD residues binding to ACE2 lie within the receptor-binding motif (RBM), which is composed by two small β-strands and 4 flexible loops. In a previous study, we pinpointed the rigidity of RBD Loop3 (L3, composed by Thr470-Pro491) as the main factor underlying the larger binding affinity of SARS-CoV-2 toward ACE2,[9] with respect to the closely related SARS-CoV. In the current set of MD simulations all the investigated systems evidence a similar RBMs flexibility, with small differences being restricted to Loop1 and 4 (L1/4, Figure S2), where N501Y is placed (Figure ).

Although not engaging direct interactions with ACE2, in the WT adduct N501@RBD intramolecularly H-bonds to Gln498, mediating the formation of a persistent H-bond network between the latter residue and Asp38@ACE2 (Figure ), thus being the most dynamically correlated residue of the whole RBM (Figure S3).[9,29] Nonetheless, in the SARBD/ACE2 and N501YRBD/ACE2 complexes, the Y501 further reinforces its prominent role in hijacking ACE2 by establishing π-stacking interactions with Tyr41@ACE2 (Figures and S4 and Table S1) and directly H-bonding to Asp38@ACE2, consistent with the N501YS-protein/ACE2 cryo-EM structure.[30] A similar H-bond pattern is also observed in μs-long MD simulations of RBD/ACE2 performed with a different force field (Table S2).[18]

Although the total binding free energies (ΔGb) of the distinct RBD/ACE2 adducts, calculated with the molecular mechanics/generalized born surface area (MM-GBSA) method,[31] do not enable one to discriminate the subtle differences between WT and mutant RBD/ACE2 adducts (Table S3), a dissection of the per-residue amino acids ΔGb contributions showed an increase, related to the N501Y substitution, with respect to the WT by 3.8 ± 2.0 and 4.4 ± 2.0 kcal/mol in SARBD/ACE2 and N501YRBD/ACE2, respectively (Figure S4), in agreement with recently reported theoretical[18] and experimental evidence.[24]

As such, N501Y, present in the highly infective UK and SA variants, possibly increases the RBD binding affinity for ACE2.[24,25] This comes along with a more effective grasping and bending of the ACE2’s α1-helix in both N501YRBD/ACE2 and SARBD/ACE2 models as compared to RBD/ACE2 (Figure ).

Figure 2

(A) Bending angle (Θ) of the angiotensin-converting enzyme 2 (ACE2)’s α1-helix defined by the Cα atoms of Phe22, Asn53, and Trp69. The receptor-binding domain (RBD), motif (RBM), and ACE2 are displayed in pink, green, and blue new-cartoons, respectively. The α1-helix@ACE2 is highlighted in silver transparent surface. (B) Distribution of Θ angle (deg) for SARS-CoV-2 and SARS-CoV RBD/ACE2[9] and (C) for SARS-CoV-2 RBD/ACE2, SARBD/ACE2, N501YRBD/ACE2, E484KRBD/ACE2, and K417NRBD/ACE2 models.

We also inspected the role of the E484K mutation common to the SA and BR variants for which no significant and reliable variation of the ΔGb could be calculated (Figures S4 and S5) in SARBD/ACE2 and E484KRBD/ACE2, respectively. Remarkably, K484 only modestly increases α1@ACE2 bending, as compared to the WT model (Figure C).

We finally assessed the role of K417,[21] whose salt-bridge with Asp30@ACE2, present in half of the RBD/ACE2 MD trajectory, is lost upon K417N mutation (Table S1). No significant variation of the ΔGb could be observed for K417NRBD/ACE2 and SARBD/ACE2 as compared to RBD/ACE2 (Figures S4 and S6),[24] and K417N does not increase the α1-helix@ACE2 bending (Figure ). Hence, the way K417N contributes to enhance the ACE2 sequestration remains elusive.

Because of the strategic location of K417, halfway of L1/4 and L3 in the RBM, tweezing the α1-helix@ACE2, we computed the cross-correlation matrices based on the Pearson’s correlation coefficient (CCs) and the per-residue sum of the cross-correlation coefficient (CCc) for the residues at the RBM/ACE2 interface (namely, we consider for each RBM residue the sum of the CCs calculated with respect to all residues of the ACE2 surface, Figure S3).[32,33] As a result, in SARBD/ACE2 and K417RBD/ACE2, the residues of the RBM exhibit the largest per-residue CCc’s. The CCs increase of SARBD/ACE2 is more marked at the L1/4 and L3 regions.

Aiming to assess whether the mutations could interfere with the RBD’s slow motions, we performed principal component (PC) analysis of WT and mutant RBD models to gather their most relevant movements (essential dynamics). As a result, PC1 and 2 of WT or all mutant RBD systems reveal the opening/closing motions of L1/4 and L3 regions, which are implicated in grasping α1-helix@ACE2 (Figure S7). Because an allosteric communication among SARS-CoV-2 mutations has been recently speculated,[34−36] we then applied dynamical network theory analysis (NWA) to decrypt the information-exchange pathways underlying the observed RBD functional dynamics and to decode whether RBD mutants can enhance the allosteric cross-talk between critical RBD’s structural elements.[32,37,38] In NWA, the protein is represented as a correlation-based weighted network. The nodes (the residues’ center of mass) are connected by edges whose numerical value (weight) indicates the correlation-strength between residue pairs (i.e., small/large weights reflect highly/poorly correlated and anticorrelated motions). By computing cross-correlations between residues along an MD trajectory, NWA finds the optimal and suboptimal signaling paths between two user-selected source (484@L3) and sink (501@L4) residues. The outcoming path lengths are thus inversely proportional to the signaling strength and to the amount of correlation existing among their tracing nodes.[37]

By performing NWA on the RBD alone we observed several cross-communication paths crossing the RBM, which in RBD/ACE2 minorly involve even K417 (Figure B,D). Remarkably, in N501Y(UK)RBD and SARBD these paths are shorter (the residues are more correlated, Figure F,H), suggesting a stronger signaling between the two RBM extremities, which may result in a more effective opening/closing of the L1/4 and L3.

Figure 3

Optimal and suboptimal signaling paths (red lines, with nodes depicted as white spheres) connecting the receptor-binding domain (RBD) residues 484 and 501, for (A) wild type (WT) RBD/angiotensin-converting enzyme 2 (ACE2), (B) WT RBD, (C) South African SARBD/ACE2, and (D) SARBD. Distribution of signaling path lengths for (E) WT RBD/ACE2 and SARBD/ACE2 adducts, (F) WT and SA isolated RBDs, (G) all investigated single mutants in the RBD/ACE2 adducts, and (H) all the isolated mutants in the RBDs.

In all RBD/ACE2 models this allosteric signaling within RBM occurs along the α1-helix@ACE2 (Figure A,C), and the path-length distribution of both SARBD/ACE2 and K417NRBD/ACE2 (Figure E,G) is shifted toward lower values. This suggests that the motions of RBM’s residues are more tightly correlated and trigger a more effective ACE2 hijacking. The similar distribution observed for SARBD/ACE2 and K417NRBD/ACE2 indicates that K417N is primarily liable for the enhanced cross-talk between critical RBD recognition loops (Figure E,G). A test of the dependence of the calculated paths from the source/sink selection (Figure S8) has been performed showing that the general trend observed in Figure E,F is maintained. However, we cannot exclude that considering a larger model of the S protein or an additional mutation present far from the RBD may perturb the observed signaling routes.[39]

To identify the residues critically involved in the signaling pathways we computed the node degeneracy (i.e., the number of times a node is present in the calculated paths). In the presence of the RBD mutations a significant variation in degeneracy is observed for those residues engaging H-bond or hydrophobic interactions at the RBD/ACE2 interface, among which are Asp38@ACE2, Asp355@ACE2, and Thr500@RBD (Figure S9 and Table S4). To further dissect the source of the increased cross-talk in the mutant RBD/ACE2 complexes we inspected whether RBD mutations alter the intra-RBD H-bonds network at the RBD/ACE2 interface. Interestingly, the main differences among the investigated systems are localized on L3–4, near the mutation sites (Table S5). In particular, a decrease of the intramolecular H-bond persistence of Asn487, which strongly H-bonds to ACE2 (Table S1), results in a higher node-degeneracy (i.e., a relevant role along the signaling route) (Table S4).

Complementarily, structural X-ray and Cryo-EM studies elucidated that K417 and E484 RBD residues establish H-bonds with distinct mAbs isolated from COVID-19 patients’ sera. Hence, K417N and E484K substitutions alter the electrostatic complementarity between the RBD and class 1 and 2 mAbs, respectively (Figure and Table S6),[40−45] impairing mAbs binding and contributing to viral escape from vaccine/disease-induced immunity.[2,41]

Figure 4

Binding mode of class 1 and class 2 (yellow and gray surfaces) monoclonal antibodies (mAbs) isolated from patients to the Spike’s receptor-binding domain (RBD), showing the use of different epitopes. The receptor-binding motif (RBM) and RBD are shown as green and pink new-cartoons, respectively. Insets disclose the key intermolecular interactions established between (from left to right) the mAbs COVA2-39,[40] C002,[41] P2B-2F6,[42] B38,[43] C105,[44] and COVA2-04.[40]

In summary, aiming to dissect the molecular basis for the higher infectivity and transmissibility of emerging SARS-CoV-2 variants, we have assessed the impact of the SA set of RBD mutations, considering them either concurrently or singularly. As a result, we disclose that while N501Y (hallmark of the UK variant) enhances the binding affinity toward ACE2 and increases the α1-helix@ACE2 bending, the SA strain exploits a two-pronged strategy to more effectively infect the host cells by (i) increasing the allosteric signaling among the pivotal RBM loops, which acting as a tweezer more effectively grasp/bend α1-helix@ACE2 and (ii) hindering the interactions with class 1 and 2 mAbs (K417N and E484K, respectively) extracted from COVID-19 patients’ sera (Figure and Table S6). Stunningly, the main actor in modulating the allosteric cross-talk among the RBD mutants appears to be K417N, whose role has remained so far elusive. In this scenario, it is tempting to argue that the BR variant, differing from the SA one only by the K417T@RBD substitution, may exploit the same strategy to foster viral propagation. Our outcomes contribute to decrypting at the atomic level the evolutionary strategies underlying the increased SARS-CoV-2 infectivity and spreading of emerging variants, setting a conceptual basis to devise next-generation therapeutic strategies against current and future viral strains.