SARS-CoV-2 B.1.617 Indian variants: Are electrostatic potential changes responsible for a higher transmission rate?

Lineage B.1.617+, also known as G/452R.V3 and now denoted by WHO with the Greek letters δ and κ, is a recently described SARS-CoV-2 variant under investigation first identified in October 2020 in India. As of May 2021, three sublineages labeled as B.1.617.1 (κ), B.1.617.2 (δ), and B.1.617.3 have been already identified, and their potential impact on the current pandemic is being studied. This variant has 13 amino acid changes, three in its spike protein, which are currently of particular concern: E484Q, L452R, and P681R. Here, we report a major effect of the mutations characterizing this lineage, represented by a marked alteration of the surface electrostatic potential (EP) of the receptor-binding domain (RBD) of the spike protein. Enhanced RBD-EP is particularly noticeable in the B.1.617.2 (δ) sublineage, which shows multiple replacements of neutral or negatively charged amino acids with positively charged amino acids. We here hypothesize that this EP change can favor the interaction between the B.1.617+ RBD and the negatively charged ACE2, thus conferring a potential increase in the virus transmission.

INTRODUCTION

The coronavirus disease 19 (COVID‐19), caused by the new coronavirus SARS‐CoV‐2, continues to spread worldwide, with more than 163 million infections and about 3.5 million deaths as of May 17, 2021 (www.who.int). To fight this dreadful disease, several safe and efficacious vaccines against SARS‐CoV‐2 are being used with remarkable effectiveness in some countries and limited availability in others. In particular, the capacity of some countries to halt SARS‐CoV‐2 spread is still limited due to inadequate resources and vaccination infrastructures. , In this scenario, several SARS‐CoV‐2 variants have been identified and have become a global threat. Some of them have been classified as variants of concern (VOCs) due to their mutations in the S1 subunit of the spike (S) protein, particularly in its receptor‐binding domain (RBD). , , One of them, identified as B.1.1.7 (α), also known as the UK variant, bears a substitution of asparagine with tyrosine on the position 501 and deletion of two amino acids in the position 69–70 of the S1 subunit. This variant has quickly spread in several European countries to become globally dominant. Other VOCs have been isolated in South Africa and Brazil and have been studied for their enhanced contagiousness and resistance to neutralization by antibodies from convalescent and vaccine‐recipient subjects , , (Table 1). Quite recently, a new variant under investigation (VUI) has been isolated from Maharashtra, India, in a setting of the highly diffusive epidemic with devastating proportions. This variant, identified as B.1.617, carries several non‐synonymous mutations. Two of them, the E484Q (or the P478K) and the L452R, are located in the RBD region, and they are critical sites for the binding with ACE2. Initial data suggest these mutations could confer increased transmission and immune evasion. , , Currently, B.1.617 comprises three subvariants, B.1.617.1‐3, with different distribution of the mutations P478K and E484Q (Public Health England). Here, we focus on biochemical and biophysical changes conferred to the B.1.617+ VUI by the P478K and E484Q mutations. We then compare these changes with other VOCs to establish whether and to what extent those amino acid changes can influence the interaction of the spike protein with ACE2, thus potentially affecting SARS‐CoV‐2 transmission and immune‐escape properties.

Table 1

Results of DynaMut and DUET analysis

Mutant sites	DynaMut						DUET
	ΔΔG (kcal/mol)a			ΔΔS (kcal/mol K)b			ΔΔG (kcal/mol)a
	6M17	6M0J	6XC4	6M17	6M0J	6XC4	6M17	6M0J	6XC4
N501Y	0.089	−0.272	0.089	−0.162	0.138	−0.312	−0.297	−0.474	−0.391
K417N	−0.399	−0.651	−1.174	0.487	0.659	0.347	−0.513	−1.295	−0.990
K417T	−0.152	−0.566	−0.832	0.198	0.507	0.193	−0.854	−1.343	−1.119
E484K	0.087	−0.101	−0.109	−0.075	0.171	0.336	0.656	0.128	0.348
E484Q	0.399	−0.644	−0.755	−0.084	0.151	0.189	0.099	−0.438	−0.319
L452R	−0.417	−0.319	−0.462	0.150	0.059	−0161	−0.548	−0.741	−0.661
T478K	−0.334	1.003	0.257	0.111	−0.385	−0.152	0.109	−0.024	0.037

Free energy difference (ΔΔG) between wild‐type and mutant structure. Negative values indicate destabilization.

Vibrational entropy difference (ΔΔS) between wild‐type and mutant structure. Positive values indicate increased flexibility.

METHODS

To perform a robust analysis three different three‐dimensional structures of SARS‐CoV‐2 spike glycoprotein have been downloaded from Protein Data Bank with the following characteristics:

SARS‐CoV‐2 RBD in complex with neutralizing antibody CC12.3 (PDB code: 6XC4, chain A),

SARS‐CoV‐2 RBD in complex with ACE2 (PDB code: 6M0J, chain E),

SARS‐CoV‐2 RBD in complex with ACE2‐B0AT1 (PDB code: 6M17, chain F).

The DynaMut server has been used to predict the impact of mutations on protein stability analyzing the folding free energy (ΔΔG) and the vibrational movement (ΔΔS), two crucial characteristics of the function and the molecular recognition of the protein. DynaMut automatically provides also results from DUET analysis, a method that combines two complementary approaches (mCSM and SDM) in a consensus prediction of ΔΔG. DUET results also have been considered in combination with DynaMut for the characterization of the mutants. PyMol (PyMol, version 2.4) suite was utilized for in‐silico mutagenesis and its Adaptive Poisson Boltzmann Solver (APBS) plugin has been used to calculate the electrostatic potential of the wild‐type and VOCs SARS‐CoV‐2 spike receptor‐binding domain (RBD) and evaluate potential differences in terms of molecular interaction with ACE2 receptor. The results have been reported within a range between −5 and +5 kT/e.

RESULTS

Protein stability

In silico prediction of the mutation impact on the RBD stability has been carried out with DynaMut. Three alternative RBD structures denoted by the PDB codes 6M17, 6M0J, and 6XC4 have been tested. These structures display small differences in the conformation of loops, especially in the one inside the receptor‐binding motif (RBM). According to the parameters of our in silico experiments, the output of DynaMut for mutant sites within loops differs depending on the starting loop conformation. For this reason, we used a normalized procedure, whereby the same mutation has been tested in each of the three different RBD structures. The effect on stability, de‐ or stabilization, has been evaluated following a majority criterion over the results of DynaMut and DUET. Detailed results have been reported in Table 1.

Starting with position 501 within a loop in the RDB interacting with the ACE2 receptor, we note that the mutation N501Y does not show any clear structural effect as the negative ΔΔG values are very close to 0.0 kcal/mol in DynaMut. At variance with DynaMut, DUET consistently assigns a destabilizing effect. For the majority rule, this mutation should be considered destabilizing (Table 1). Also, position 417 is within the interface α‐helix where Lys interacts with ACE2 Asp30. However, in this case, our results predict that both mutants (K417N and K417T) could destabilize the protein and increase local flexibility. Similarly, Glu484 is in an interfacial loop and interacts with ACE2 Lys31. Our data predict that mutation E484K is stabilizing, although data from DynaMut and DUET do not match. E484Q exerts a destabilizing effect, again with no strong consensus from the two methods. We note that glutamic acid is a polar, negatively charged, hydrophilic amino acid while lysine is a basic, charged and partly aliphatic amino acid, and its ε‐amino group often participates in hydrogen bonding, salt bridges, and covalent interactions. Both mutants E484K and E484Q would likely increase local molecular flexibility. L452 is in a short β‐strand, and it is exposed to the solvent. Apparently, it does not interact directly with ACE2. Mutation L452R is predicted to be destabilizing with increased local flexibility. Position 478 is in a loop in the proximity of the ACE2 although not in direct contact with it. Mutation T478K is stabilizing and is predicted to decrease local protein flexibility.

Surface and interface analysis

The entire set of mutations in positions 417, 452, 478, 484, and 501 are in the spike RBD at the interface with the ACE2 N‐terminal helix (Figure 1). However, we note that the mutant positions 452, 478, 484, and 501 are within the RBM, containing residues that bind to ACE2. In contrast, mutant position 417 is located outside the motif. According to our analysis, the mutations in positions 417, 484, and 501 might increase the spike binding affinity with the ACE2 receptor. In particular, the Tyr replacing Asn501 may form an aromatic interaction with ACE2 Tyr41, a hydrogen bond with ACE2 D38, and a potential cation–π interaction with ACE2 Lys353. In addition, substitutions of Glu484 with Lys or Gln may form hydrophobic interactions to Ile472, Gly482, Phe486, Cys488, and Tyr489. Our data also indicate that replacing Lys or Gln with Glu484 abolishes the interfacial salt bridge between Glu484 and ACE2 Lys31. Due to the fact that Lys417 is solvent‐exposed and forms salt‐bridge interactions with Asp30 of ACE2, the replacement of Thr/Asn with Lys417 could abolish this interaction. Moreover, although the mutations in positions 452 and 478 are within the receptor‐binding motif, our analysis does not show a direct interaction with ACE2.

Figure 1

Comparison between the wild‐type and variant spike receptor‐binding domains (RBDs). Protein surface is colored according to the electrostatic potential. Color scale ranges from −5 kT/e (red) to +5 kT/e (blue) as reported by the bar under the wild‐type RBD. Position of the mutant sites is indicated by a circle and an attached label. Red arrows mark the area of increased positive potential in the RBD Indian variants

Electrostatic potential

We note that a major, global effect of the mutations characterizing the Indian variants is represented by the alteration of the RBD surface electrostatic potential. In particular, in the lineage B.1.617.1 (κ) the uncharged and hydrophobic residue Leu452 changes to the positively charged residue Arg, and the negatively charged residue Glu484 is replaced by the uncharged residue Gln. In contrast, the B.1.617.2 (δ) lineage shares the same mutation in position 452, but it has another mutation in position 478 where the neutral residue Thr changes to the positively charged Lys. The presence of two positively charged residues in the variant B.1.617.2 (δ) increases the positive electrostatic potential surface (Figure 2).

Figure 2

Ribbon model of the interface between ACE2 (deep teal) and receptor‐binding domain (RBD) (orange). Sidechains of relevant residues are displayed as stick models and labeled. The two mutations at the RBD sites 417 and 484 have been simultaneously displayed

DISCUSSION

Given that SARS‐CoV‐2 variants are characterized by evolutionary and genetic changes accumulated in the genome, the use of new and improved phylodynamic techniques for the study of how epidemiological, immunological, and evolutionary processes contribute to shaping viral phylogenies is extremely important and useful. Here, using a well‐structured software workflow, we provide evidence of a strong system able to carry out a quick, systematic, and reproducible screening of the SARS‐CoV‐2 genome isolates. Protein mutations were identified by scanning high‐quality SARS‐CoV‐2 genomes variants downloaded from the GISAID databank. We reasoned that the most likely health‐threatening properties of SARS‐CoV‐2 VOC rely on fine biochemical and biophysical changes that eventually impact the RBD interaction with the ACE2 receptor present on the host's cell surface. We thus characterized B.1.617+ SARS‐CoV‐2 VUI using in‐silico methods capable of predicting the effect of mutations on S‐RBD protein stability and its electrostatic potential. The B.1.617+ mutations have been investigated by comparing most of the already known mutations of previously reported VOCs. To further justify our assumption, we note that the intensely investigated D614G substitution of the spike protein, early reported in Italian isolates, , , and subsequently attributed with increased virus transmissibility was found to enhance the protein torsional ability and potentially affecting its stability. Regarding the commonly denominated Indian variants, these constitute the new SARS‐CoV‐2 lineage B.1.617, which is actually composed of a family of three subvariants, namely B.1.617.1 (κ), B.1.617.2 (δ), and B.1.617.3. This lineage, which emerged in India in October 2020, has since spread to other countries, particularly in the United Kingdom, in settings with a high density of resident immigrants from India. Data from Public Health England registry shows that the subvariant B.1.617.2 (δ) has become epidemiologically prevalent.  Due to its subvariants composition, the variant has been designed as B.1.617+ by the WHO.

A dichotomic behavior is observed with variant B.1.617.2 (δ). In contrast, it lacks mutation E484Q which is present in the other two lineages and was initially suspected to confer a certain degree of resistance to antibody neutralization. In contrast, this subvariant has a mutation at site 478 where a lysine replaces the proline. Of note, variant B.1.617.2 (δ) is indeed characterized by a major shift toward increased positive electrostatic potential because of three amino acid changes from negative or neutral to clearly positive charge, as shown in Section 3.

Subvariants B.1.617.1 (κ) and B.1.617.3 have both the double mutations E484Q and the L452R. Although initially believed to enhance the antibody‐escape potential, it has been shown that the B.1.617.1 (κ) subvariant is pretty neutralized by the majority of sera from convalescent individuals and all sera from vaccinated subjects. Nonetheless, this variant has been shown to be more pathogenic than the B.1.1.7 (α) variant in an experimental model of SARS‐CoV‐2 infection in hamsters.

Our data indicate that most of the mutants are predicted to destabilize the RBD structure, except for E484K and T478K. It is conceivable that destabilization alters binding affinity to ACE2 and to neutralizing antibodies. As noticed above, the influence of the mutations on the spike surface properties is evident in the B.1.617+ lineage, particularly in the B.1.617.2 (δ) lineage subvariant, where the positive electrostatic surface potential is markedly increased. This may favor RBD interaction with the negatively charged ACE2, which, in turn, would then increase affinity for the ACE2 receptor. All of these changes have the potential to eventually modify infectivity, pathogenicity, and virus spread. Regarding differential binding to neutralizing antibodies, previous studies suggested that VOCs RBD changes in the electrostatic potential surface could induce SARS‐CoV‐2 antibody evasion, and even single amino acid changes that marginally affect ACE2 affinity could greatly influence nAbs affinity. Several factors have been demonstrated to affect the impact of VOCs. For example, it has been observed an increased effect at pH associated with nasal secretions (from 5.5 to 6.5). For this reason, additional experiments both in vitro and in vivo are needed to establish the biological significance of SARS‐CoV‐2 mutations and how the interactions between mutations and local cellular microenvironment influence the clinical outcome and the transmission dynamics of this virus.

AUTHOR CONTRIBUTIONS

Stefano Pascarella and Martina Bianchi performed the in‐silico analysis and proposed data interpretation. Antonio Cassone and Massimo Ciccozzi proposed addressing the issue of this study, wrote the initial draft of the manuscript, and contributed to writing the final version of the manuscript. Stefano Pascarella, Martina Bianchi, Antonio Cassone, Massimo Ciccozzi, Davide Zella, Francesca Benedetti, Francesco Broccolo, Roberto Cauda, Arnaldo Caruso, Silvia Angeletti, Marta Giovanetti, and Domenico Benvenuto discussed the data, their interpretation, and contributed to writing the final version of the manuscript. All Authors have seen and approved the initial draft and the revised, final text of the manuscript.