Characterization of mutations modulating enhanced transmissibility of SARS-CoV-2 B.1.617+ (Delta) variant using In Silico tools.

Since the beginning of the of SARS-CoV-2 (Covid-19) pandemic, variants of concern (VOC) have emerged taxing health systems worldwide. In October 2020, a new variant of SARS-CoV-2 (B.1.617+/Delta variant) emerged in India, triggering a deadly wave of Covid-19. Epidemiological data strongly suggests that B.1.617+ is more transmissible and previous reports have revealed that B.1.617+ has numerous mutations compared to wild type (WT), including several changes in the spike protein (SP). The main goal of this study was to use In Silico (computer simulation) techniques to examine mutations in the SP, specifically L452R and E484Q (part of the receptor binding domain (RBD) for human angiotensin-converting enzyme 2 (hACE2)) and P681R (upstream of the Furin cleavage motif), for effects in modulating the transmissibility of the B.1.617+ variant. Using computational models, the binding free energy (BFE) and H-bond lengths were calculated for SP-hACE2 and SP-Furin complexes. Comparison of the SP-hACE2 complex in the WT and B.1.617+ revealed both complexes have identical receptor-binding modes but the total BFE of B.1.617+ binding was more favorable for complex formation than WT, suggesting L452R and E484Q have a moderate impact on binding affinity. In contrast, the SP-Furin complex of B.1.617+ substantially lowered the BFE and revealed changes in molecular interactions compared to the WT complex, implying stronger complex formation between the variant and Furin. This study provides an insight into mutations that modulate transmissibility of the B.1.617+ variant, specifically the P681R mutation which appears to enhance transmissibility of the B.1.617+ variant by rendering it more receptive to Furin. Crown

Introduction

In December 2019, a novel coronavirus, SARS-CoV-2 (referred to as COVID-19), was reported in Wuhan, China. Since then, several variants of concern (VOC) have emerged that can be more virulent, transmissible and increase the risk of death (Islam et al., 2022). Furthermore, variants have potential to decrease effectiveness of current mitigation strategies and risk public health globally (Babuadze et al., 2022). VOCs emerge due to mutations in key genes needed for pathogenesis. As a result, mutations have a crucial role in the ongoing evolution and emergence of novel SARS-CoV-2 variants (Shahhosseini et al., 2021a). As such, a new variant of SARS-CoV-2 was detected in India in October 2020, and is the cause of a second deadly wave of the COVID-19 pandemic in this country. Gene mapping of SARS-CoV-2 B.1.617+ (Delta/Indian variant) variant revealed several amino acid mutations in both the nucleocapsid and the spike protein (SP) (Gupta, 2021). The SP had 8 amino acid mutations in B.1.617+, specifically, G142D, E154K, L452R, E484Q, D614G, P681R, Q1071H, H1101D. Among the amino acid mutations in the SP, L452R and E484Q occur in the receptor binding domain (RBD), and P681R occurs in the protease cleavage site (Gupta, 2021).

The SP of SARS-CoV-2, which is critical for transmissibility of the virus, is comprised of subunit S1 (amino acids from 14 to 685) and subunit S2 (amino acids from 686 to 1147). The RBD is part of the S1 subunit and recognizes human angiotensin-converting enzyme 2 (hACE2) on the surface of host cells and hence is responsible for virus attachment. As a result, any mutation present in the RBD may potentially affect binding affinity to hACE2, increasing or decreasing SARS-CoV-2 transmissibility. The S2 subunit contains protease cleavage sites at positions 685, 696, and 815, which can be recognized by proteases (e.g. Furin, Transmembrane protease serine 2, Cathepsin L, and Trypsin). A vital step in SP activation and virus fusion into the host cell is cleavage of the virus by proteases. It was previously shown that a unique characteristic of SARS-CoV-2 is the presence of four residues (bold) that are part of the Furin recognition cleavage motif (underline) ‘681 P /SV687’. Furthermore, previous studies demonstrated that although the multi-basic cleavage motif ‘RRAR’ is essential for efficient cleavage of the SP by Furin, six residues upstream and two residues downstream of the cleavage site have a significant role in modulating Furin cleavage efficiency (Shiryaev et al., 2013; Örd et al., 2020). The ability of Furin to recognize the Furin cleavage sites at the S1/S2 boundary is essential for virion fusion and entry into human cells and any mutations in the protease cleavage site may have an influence SARS-CoV-2 transmissibility (Wu et al., 2020).

Two methods can be used to examine the effect of specific mutations using In Silico models is to determine the binding free energy (BFE) and the length of the hydrogen bonds between complexes. The contribution of protein residues in complex formation with a ligand is represented by the BFE. As a result, the protein-ligand binding affinities can be predicted using BFE calculations. The strength of the bonds between the virus residues and the receptors on the surface of host cell can potentially be used to assess virus transmissibility due to increased binding affinity between a virus protein and the host cell receptors (Zou et al., 2020; Maurya et al., 2020). Moreover, the length of the hydrogen bond (H-bond) between viral and host amino acids has been shown to affect the probability of molecular interactions, which are essential for the formation of a stable complex and virus attachment to the host cells (Jakhmola et al., 2021). The bond length of 2.2–2.5 Angstrom (Å), 2.5–3.2 Å, and 3.2–4.0 Å respectively are the cut-offs for a strong, moderate, and weak H-bond (Jeffrey and Jeffrey, 1997; Marqusee and Baldwin, 1987).

The main goal of this study was to use In Silico approaches to characterize three mutations of the B.1.617+ variant of SARS-CoV-2 present in the SP, specifically, L452R, E484Q and P681, and investigate the role of these mutations in the emergence of a more contagious variant.

Materials and methods

Split tree of SARS-CoV-2 B.1.617+ variant

A total number of 95 full-length SARS-CoV-2 sequences representing all clades were retrieved from the global initiative on sharing all influenza data (GISAID) database (https://www.gisaid.org/), similar sequences from the same origin and clade were excluded from final analysis. Additionally, the first SARS-CoV-2 sequence from Wuhan was retrieved from GenBank (https://www.ncbi.nlm.gov) as the WT strain (reference sequence: NC_045512). The split network of SARS-CoV-2 variants was constructed by EqualAngel method using the SplitsTree 4.0 software (Huson and Bryant, 2006; Shahhosseini et al., 2016; Shahhosseini et al., 2021b).

Modeling In Silico docking of WT or B.1.617+ variant SP with hACE2/Furin

The 3D structures for each protein including the WT SP (NCBI accession: NC_045512), the B.1.617+ variant SP (GISAID: EPI_ISL_1704384), hACE2 (NCBI accession: 6m17_B), and Furin (NCBI accession: NP_002560) were constructed using the RaptorX server (http://raptorx.uchicago.edu). The accuracy of structures were re-evaluated by SWISS-MODEL workspace (http://swissmodel.expasy.org) (Schwede et al., 2003; Waterhouse et al., 2018). Subsequently, the 3D structures with the lowest root mean square deviation (RMSD) and qualitative model energy analysis (QMEAN) Z-scores around zero (indication of models with high quality) were selected to produce a Protein data bank (PDB) file for further analysis (Källberg et al., 2012; Joosten et al., 2010). High Ambiguity Driven protein-protein DOCKing (HADDOCK) v.2.2 (https://www.bonvinlab.org/software/haddock2.2/haddock-start/) was used to build interactive docking complexes between hACE2 and Furin with WT and the B.1.617+ variant SP, and the complexes were visualized by 3DBIONOTES web server (https://3dbionotes.cnb.csic.es) (Van Zundert et al., 2016; Wassenaar et al., 2012; Segura et al., 2019). Additionally, the precision of the docking mode and interactions in the built docked structures were validated using two independent servers; Fast Rotational DOCKing (FRODOCK) v.3.12 (http://frodock.chaconlab.org) and HDOCK (http://hdock.phys.hust.edu.cn) (Ramírez-Aportela et al., 2016; Yan et al., 2020). Among the tentative clusters, the docked complexes with the highest HADDOCK score, a larger reproducible cluster size, and the lowest possible RMSD, suggesting the high likelihood for accurate structures were selected for further analysis (Källberg et al., 2012; Yan et al., 2020).

BFE of WT or B.1.617+ variant SP-hACE2/Furin complex

Docked complex PDB files were used to assess the impact of mutations on the binding affinity of the SP for hACE2. The BFE between WT and B.1.617+ SP of SARS-CoV-2 and hACE2/Furin in protein-ligand docking complexes was measured using the Molecular Mechanics with Generalized Born and Surface Area (MM/GBSA) method implemented in the HawkDock web server (http://cadd.zju.edu.cn/hawkdock/) (Hou et al., 2011; Wang et al., 2019). This server uses the ATTRACT docking algorithm and the HawkRank scoring function (Weng et al., 2019).

Molecular interaction between WT or B.1.617+ variant SP with hACE2/Furin

Docked complex PDB files were used to infer protein-protein interaction (PPI) in the docked structures of hACE2/Furin with WT or the B.1.617+ variant SP using interacting residue plots created by the Dimensional reduction plot (DIMPLOT) tool implemented in Ligplot+ v.2.2.4 (EMBL-EBI, Cambridge, UK) (Wallace et al., 1995). Runtime parameters relevant to the HBPLUS program were set to 4 Å to determine possible H-bonds and non-bonded contacts in order to generate Ligplot+ diagrams. Only interactions with strong H-bonds (length of 2.5 Å or less (Jeffrey and Jeffrey, 1997)) were considered for analysis.

Prediction of Furin cleavage efficiency

The impact of mutations on Furin cleavage efficiency at the S1/S2 boundary of WT or the B.1.617+ variant SP was evaluated. For this purpose, the ProP 1.0 server (http://www.cbs.dtu.dk) was employed to compute the Furin cleavage efficiency of the multi-basic cleavage site. The minimum score for Furin cleavage was 0.5, and the maximum score was 1 (Duckert et al., 2004).

Results

Split tree of B.1.617+ variant SARS-CoV-2

Analysis with the Split network revealed that sequences of SARS-CoV-2 clustered into eight clades (Fig. 1 ). The B.1.617+ variant clustered in a distinct clade (B.1.617+) compared to previously identified VOC; clade GR including the B.1.1.7 and P.1 variants, and clade GH including B.1.351 and B.1.427. Furthermore, the B.1.617+ variant diverged into three sub-clades as B.1.617.1, B.1.617.2, and B.1.617.3 within the main B.1.617+ clade.

Fig. 1

The Neighbor-Net network of SARS-CoV-2 was constructed using the SplitsTree 4.0 software. Each clade is identified by different colors (upper left). In addition to the B.1.617+ variant, four VOC are shown, including two variants in clade GR (B.1.1.7 and P.1) and two variants in clade GH (B.1.351 and B.1.427).

Gene mapping of the three sub-clades of B.1.617+ revealed that B.1.617.1 and B.1.617.3 share two mutations in the RBD (L452R and E484Q) and one mutation in the Furin cleavage loop (P681R), whereas, B.1.617.2 lacks the E484Q mutation in the RBD and only shares the L452R and P681R mutations with the two other sub-clades.

BFE of SP-hACE2 docked complexes

Two of the mutated residues 452 and 484, in both WT or the B.1.617+ variant SP, do not appear to significantly interact with hACE2. However, as shown by plots, the E484Q mutation seems to slightly increase the residue accessibility (Fig. 2 A and D).

Fig. 2

Structures of hACE2 docked with WT (B) or the B.1.617+ variant (E) SARS-CoV-2 SP are modeled using SWISS-MODEL. The interacting residues/residue accessibility (WT; A, the B.1.617+ variant; D), and the enlargement of SP-hACE2 docked complexes (WT; C, the B.1.617+ variant; F) are computed using HADDOCK v.2.2. The contribution of mutant residues to the total BFE (kcal/mol) in the SP-hACE2 docked complexes for WT and the B.1.617+ variant are computed by the MM/GBSA method, implemented in the HawkDock web server, while mutations are shown with red boxes. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

In Silico analysis demonstrated that L452 contribution to the total BFE in the WT variant is +0.01 kcal/mol, while the R452 mutation decreased the BFE to −1.02 kcal/mol in the B.1.617+ variant, which favors complex formation. Therefore, considering BFE, it seems that the L452R mutation favors complex formation between the RBD and hACE2. Furthermore, residue E484 in the flexible loop region of the RBD in the WT which contributes to a BFE of +1.91 kcal/mol, an undesirable contribution to the total binding energy. However, E484Q mutation decreases the BFE to +0.39 kcal/mol for the B.1.617+ variant, which is not a significant contribution. The total BFE of SP-hACE2 docked complexes in SARS-CoV-2 WT and the B.1.617+ variant are −64.88 and −65.67 kcal/mol, respectively (Fig. 2).

Molecular interactions between hACE2 with WT or B.1.617+ variant SP

A total of 14 residues in both WT and the B.1.617+ variant SP are considered to be in close contact with hACE2 based on a cutoff distance of 4 Å. Therefore, amino acids at positions 417, 496, 449, 455, 456, 486, 487, 489, 493, 498, 500, 475, 502, and 505 could be contributing to the SP-hACE2 interaction.

Both residues L452 and E484 based on molecular interactions do not appear to contribute to SP-hACE2 interaction since they are located in regions with a lower probability of contact with hACE2. Results revealed that residue L452 (Leu452) in WT and R452 (Arg452) in the B.1.617+ variant have a mean distance of 9.98 Å from hACE2 (distances of 4 Å or longer in length are not shown in Fig. 3 ). Similarly, residue E484 (Glu484) in WT and Q484 (Gln484) in the B.1.617+ variant also do not appear to interact with hACE2, as the H-bond length is longer than 4 Å. Furthermore, as shown in Fig. 3, the H-bonds and hydrophobic interactions between hACE2 and WT or the B.1.617+ variant RBD have overlap with respect to docking predictions, suggesting that substantial conformational changes may not be present in RBD due to mutations in the B.1.617+ variant SP (Fig. 3).

Fig. 3

Interacting residues between WT or the B.1.617+ variant SP and hACE-2 are merged, while equivalent side chains are overlaid and marked by red circles. The DIMPLOT tool implemented in Ligplot+ v.2.2.4 are used to infer H-bond (green dots) and hydrophobic (red dots) interactions. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

BFE of SP-Furin docked complexes and cleavage efficiency of SP

The protein-protein docking result showed that residues T678-Q690 in WT and residues Q677-S691 in the B.1.617+ variant interact with Furin. Furthermore, it seems that the mutation P681R can increase the residue accessibility at the Furin cleavage loop in the B.1.617+ variant (Fig. 4 A and D; interacting residues are colored as yellow). The MM/GBSA analysis shows that the residue P681 in WT contributed to the total BFE by −3.77 kcal/mol, while the corresponding residue in the B.1.617+ variant (R681) lowered the contribution to −8.12 kcal/mol, which favors complex formation. In general, the total BFE of the multi-basic cleavage site-Furin docked complexes in WT SARS-CoV-2 is −60.32 kcal/mol, while in the B.1.617+ variant is −71.36 kcal/mol (Fig. 4).

Fig. 4

Structures of Furin docked with WT (B) or the B.1.617+ variant (E) SARS-CoV-2 SP are modeled using SWISS-MODEL. The amino acid annotation shows the interacting residues and residue accessibility (peaks indicate higher accessibility) of WT (A) or the B.1.617+ variant (D) SARS-CoV-2 for Furin. Sections A and D, as well as the enlargement of SP-Furin docked complexes (WT; C, the B.1.617+ variant; F) are computed using HADDOCK v.2.2. The BFE values (kcal/mol) are calculated using the MM/GBSA method implemented in the HawkDock web server. Mutations are indicated with red boxes. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Furthermore, the cleavage efficiency analysis using ProP demonstrated that P681R mutation increased the Furin score from 0.620 in WT (NSPRRAR/SV) to 0.698 in the B.1.617+ variant (NSRRRAR/SV), suggesting increased cleavage efficiency of Furin (Fig. 5 ).

Fig. 5

The Furin cleavage efficiency of WT (A) and the B.1.617+ variant (B) were predicted by ProP 1.0 server.

Molecular interactions between Furin with WT or B.1.617+ variant SP

The active site pocket of Furin includes residues R185, M189, D191, N192, R193, E229, V231, D233, D259, K261, R298, W328, and Q346, and has a catalytic role. The binding pocket of Furin accommodates the Furin cleavage loop of SARS-CoV-2 SP, and it can strongly interact with residues 678TNSPRRAR/SVASQ690 of WT SP and 677QTNSRRRAR/SVASQS691 in the B.1.617+ variant SP (H-bonds of 4 Å or longer in length are not shown in Fig. 6 ).

Fig. 6

Interacting residues between WT (A) or the B.1.617+ variant (B) SP and Furin are shown, while red circles mark equivalent side chains. The DIMPLOT tool implemented in Ligplot+ v.2.2.4 was used to infer H-bond (green dots) and hydrophobic (red dots) interactions. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

As shown in Fig. 6, the H-bonds and hydrophobic interactions between Furin and WT or the B.1.617+ variant SP are identical, except at residue 681. In the WT, the residue P681 (Pro681) forms a hydrophobic interaction with G229 (Gly229) and V231 (Val231) on Furin, whereas the residue R681 (Arg681) in the B.1.617+ variant forms an H-bond interaction with residue G229 (Gly229) on Furin, which is a strong bond with the length of 2.81 Å (Fig. 6).

Discussion

The B.1.617+ variant with mutations present at L452R and E484Q in the RBD and P681R in the Furin cleavage site of the SP is thought to be more transmissible according to available epidemiological data (Hart et al., 2022). Although the exact dissemination of B.1.617+ is uncertain, 53.28 % (2325 from 4363) of all samples submitted to the GISAID from B.1.617+ with collection dates between March 1st (starting the second wave of pandemic) and May 22nd (data set preparation for this study) in 2021 are classified as B.1.617+. According to the Split tree analysis and mutation signature, the B.1.617+ variant has emerged separately in India. Split Network demonstrated that the B.1.617+ variant diverged into three sub-clades. B.1.617.2 corresponds to 74.37 % (7504 from 10,090) of all the B.1.617+ variants submitted to GISAID followed by B.1.617.1 (24.76 %) and B.1.617.3 (0.86 %). Subsequently, this study aimed to predict the binding affinity of the SP-hACE2 and the SP-Furin complex in these emerging variants from India (B.1.617+).

Overall, examining the SP-hACE2 complex binding energy determined that neither L452R or E484Q mutations contribute significantly to the total binding energy but slight differences between the WT and B.1.617+ are present. Specifically, the BFE calculation demonstrated that residue L452 in the WT variant and R452 in the B.1.617+ variant contribute to the total binding energy by +0.01 kcal/mol and −1.02 kcal/mol, respectively. The values moderately favor complex formation between the SP and hACE2 in the B.1.617+ variant by lowering the binding energy. Furthermore, the mutation E484Q has the potential to transform the BFE from an undesirable contribution in the WT variant (+1.91 kcal/mol) to a nearly neutral contribution in the B.1.617+ variant (+0.39 kcal/mol) (Fig. 2).

The BFE results alone could not be used to predict the binding affinity of the SP for hACE2, therefore, the feasibility of interaction between the RBD of the SP and hACE2 was evaluated by measuring the length of H-bonds at the interface sites. Molecular interaction between hACE2 and WT or the B.1.617+ variant RBD showed that neither L452R nor E484Q were found within the key contact residues of the RBD, which is consistent with a prior study (Celik et al., 2021). Thus, residues 452 or 484 did not appear to form stable H-bond with hACE2 (bond length of <4 Å). Moreover, the L452R and E484Q mutations did not cause significant structural changes in the RBD (for the key contact residues) as demonstrated by the molecular interaction models (Fig. 3). Moreover, these mutations (L452R and E484Q) do not have a significant impact on the SP-hACE2 binding affinity and most likely do not modulate the hypothesized enhanced transmissibility of the B.1.617+ variant. The results are consistent with epidemiological evidence, as B.1.617.2, the most common of the three B.1.617+ sub-clades, lacks the E484Q mutation providing evidence that this mutation is not needed for enhanced transmission of the variant. This finding is in agreement with another study, which found that the binding affinity of the B.1.617+ variant RBD to hACE2 did not show any substantial increase due to mutations in RBD, suggesting that there must be other alternative mechanisms or mutations causing enhanced transmissibility of the B.1.617+ variant (Kumar et al., 2021). Apart from binding affinity, it is worth noting that amino acid substitutions at positions L452R and E484Q could have a role in allowing the virus to escape neutralizing antibodies (Villoutreix et al., 2021; Starr et al., 2021).

Determining that L452R and E484Q does not appear to substantially increase the binding affinity of the B.1.617+ variant to hACE2, it was hypothesized that the increased transmissibility of the B.1.617+ variant (based on epidemiological observations) is likely due to the P681R mutation in the upstream of the Furin cleavage site (R685/S686). Thus, in order to assess other mechanisms that can elucidate increased B.1.617+ transmissibility, the Furin cleavage loop at the S1/S2 boundary was explored. Study results demonstrated, firstly, the total BFE was reduced from −60.32 kcal/mol in WT to −71.36 kcal/mol in the B.1.617+ variant, implying that the P681R mutation in the cleavage loop renders it more accessible to Furin, thus, enhancing the infection efficiency (Fig. 4). Secondly, the predicted cleavage efficiency (Furin score) at the S1/S2 boundary of the B.1.617+ variant (0.698) is higher than for WT (0.620) (Fig. 5), suggesting that the SP is more easily cleaved by Furin, further facilitating virus fusion with host cell membranes. Thirdly, the residue P681 in the WT virus makes a hydrophobic interaction with G229 and V231 on Furin, whereas the residue R681 in the B.1.617+ variant makes a strong H-bond interaction with residue G229 on Furin (bond length of 2.81 Å) (Fig. 6). Altogether, the results from this study indicate that the enhanced transmissibility of the B.1.617+ variant could be due to the P681R mutation at the upstream of the Furin cleavage site.

In Silico investigations from this study are in accordance with a previous study which demonstrated that the addition of an R (Arg) to the position 681 jointly with K (Lys) at position 684, increased cell-cell spread of the SARS-CoV-2 and potentially changed the pathogenicity (Hoffmann et al., 2020). Similarly, the amino acid sequences ‘RRRKR/GL’ has also been identified in the Furin cleavage site of H5N8 (Örd et al., 2020).

Apart from the computational data presented here, which introduces P681R as a modulating mutation for the potentially enhanced transmissibility of the B.1.617+ variant, epidemiological data from Uganda identified the emergence of a new variant (Clade S, lineage A.23.1) with the P681R mutation. This new variant has dominated the Uganda pandemic, implying the likelihood of augmented transmission (Bian et al., 2021). Combining the modeling data from the current study with epidemiological data from Uganda and India support the hypothesis that the P681R mutation is linked to augmented transmissibility.

The UK variant, B.1.1.7, has a P681H mutation in the Furin cleavage loop, which is thought to play a role in the enhanced transmissibility of this variant (Davies et al., 2021; Rambaut et al., 2020). The cleavage efficiency analysis for the UK variant with the P681H mutation in the cleavage loop (NSHRRAR/SV) was 0.701 (Supplementary Fig. 1), implying that the UK variant has the most efficient Furin cleavage of the virus, followed by the B.1.617+ variant with R681 (0.698) and WT with P681 (0.62) (Fig. 5). This finding that mutations in the Furin cleavage loops of both the UK (B.1.1.7) and B.1.617+ support more efficient cleavage, combined with same type of mutation (as B.1.617+) noted in the dominant Uganda strain, strongly supports the theory that mutations in the Furin cleavage loop have an important role in variant transmissibility.

Supplementary Fig. 1

The Furin cleavage efficiency of B.1.1.7 (UK) was predicted by ProP 1.0 server.

Although, the main focus of this study was to study mutations in the SARS-CoV-2 SP and their potential impact on altering the binding affinity of the B.1.617+ variant to hACE2 and Furin, it should be noted that virus transmissibility is a multifactorial mechanism. Biological factors such as the state of the host immune system (which varies by age and/or diagnosis of other diseases) and epigenetics (which affects polymorphism in human genes encoding SARS-CoV-2 receptors (e.g. hACE2 or TMPRSS2) or cleavage proteases (e.g. Furin) on the surface of human cells) are significant determinants in susceptibility to a viral variant. These factors may alter binding affinity for the viral SP, and thereby enhancing/weakening the transmissibility of a viral variant in human populations (Hou et al., 2020). Furthermore, it is worth mentioning that non-biological factors like human social activity (e.g. avoiding event gatherings, respecting social distancing and lockdowns) and/or utilizing personal protective measures (e.g. wearing a mask and using sanitizer) can all play a role in preventing the rapid spread of a variant in a given time and space.

Taking genetic and epidemiological data together, VOC can be defined when the augmented transmissibility and deteriorating epidemiological situations are associated with mutation signatures known to enhance the binding affinity of the virus (e.g. N501Y) (Shahhosseini et al., 2021c). Accordingly, the most contagious VOC until now are lineages P.1 (Brazil) (Coutinho et al., 2021), B.1.1.7 (UK) (Volz et al., 2021; Washington et al., 2021) and B.1.351 (South Africa) and since the majority of SARS-CoV-2 cases detected in India's second wave of the pandemic were the B.1.617+ variant (epidemiological data) and this variant contains major mutation signatures such as the P681R mutation (genetic data), this study would consider the B.1.617+ variant as a VOC. Further experimental studies are required to confirm this statement.

Investigating the modulating factors of Furin cleavage site and RBD binding affinity to hACE2 paves the way towards developing novel therapies that inhibit the Furin cleavage site or RBD of SARS-CoV-2, along with predicting the transmissibility of a variant using computational modeling at the protein level.

Conclusion

Since the beginning of the SARS-CoV-2 pandemic, mutations have played a crucial role in the ongoing evolution and emergence of novel SARS-CoV-2 VOC leading to several deadly waves of infection. These mutations often increased the transmissibility of SARS-CoV-2 and challenged healthcare systems worldwide. The results from this study focused on three different mutations on the SP of the B.1.617+. Overall, the results introduced compelling evidence that the enhanced transmissibility of the B.1.617+ SARS-CoV-2 variant is mainly due to P681R mutation on the upstream cleavage motif rather than the L452R or E484Q mutations present in the RBD. Future studies aimed at examining other identified mutations may further elucidate the transmissibility of B.1.617+ variant using computational modeling. In Silico models have the advantage of analyzing multiple data sets quickly without the need for live virus and can assist in guiding future in vitro/in vivo experimentation by focusing research on areas most likely to be responsible for the observed effect and leading to enhanced disease mitigation.

The following are the supplementary data related to this article.

Supplementary Table 1

We gratefully acknowledge the following Authors from the Originating laboratories responsible for obtaining the specimens and the Submitting laboratories where genetic sequence data were generated and shared via the GISAID Initiative, on which this research is based.

Funding

This study was not supported by a specific fund.

CRediT authorship contribution statement

Nariman Shahhosseini: Conceived and designed the study, Data collection, Analysed the data, Wrote the manuscript.

Declaration of competing interest

No conflict of interests to declare.