Identification of novel mutations in RNA-dependent RNA polymerases of SARS-CoV-2 and their implications on its protein structure

Sequence retrieval

We downloaded all SARS-CoV2 sequences from the NCBI virus database as shown in Table 1. As of now, there are 28 SARS-CoV-2 sequences from India have been deposited in this database, out of which, two sequences are not complete; therefore, we retrieved 26 datasets of SARS-CoV-2 (samples are from Indian). As a reference, we downloaded the sequence of SARS-CoV-2 that was first reported genome sequence deposited in the NCBI virus database from the ‘Wuhan wet sea food market area’ from the early days of COVID-19 pandemic (Wu et al., 2020). This virus was formerly called ‘Wuhan seafood market pneumonia virus’ with the accession number YP_009724389.

Sequence alignments and structure

All the RdRp protein sequences were aligned by multiple sequence alignment platform of CLUSTAL Omega (Madeira et al., 2019). Clustal Omega is a multiple sequence alignment programme that uses seeded guide trees and HMM profile-profile techniques to generate alignments between three or more sequences. The alignment file was carefully studied and differences in the amino acid changes were recorded.

Secondary structure predictions

We used CFSSP (Ashok Kumar, 2013) (Chou and Fasman secondary structure prediction), an online server, to predict secondary structures of SARS-CoV-2 RdRp protein. This server predicts the possibility of secondary structure such as alpha helix, beta sheet, and turns from the amino acid sequence. CFSSP uses Chou-Fasman algorithm, which is based on analyses of the relative frequencies of each amino acid in secondary structures of proteins solved with X-ray crystallography.

RdRp dynamics study

To investigate the effect of mutation on the RdRp protein structural conformation, its molecular stability and flexibility, we used DynaMut software (Rodrigues, Pires & Ascher, 2018) (University of Melbourne, Australia). To run this software we first downloaded the known protein structure of SARS-CoV-2 RdRp from RCSB (PDB ID: 6M71 (Gao et al., 2020), 7BV1 and 7BV2 (Yin et al., 2020)) and used it for analysis. Next, the 6M71, 7BV1 and 7BV2 structure was uploaded on DynaMut software and effect of mutation in various protein structure stability parameters such as vibrational entropy; the atomic fluctuations and deformation energies were determined. DynaMut, a web server implements well established normal mode approaches that is used to analyse and visualise protein dynamics. This software samples conformations and measures the impact of mutations on protein dynamics and stability resulting from vibrational entropy changes and also predicts the impact of a mutation on protein stability.

Identification of mutations in RdRp protein present in Indian isolates

The SARS-COV-2 sequencing data was downloaded from NCBI (NCBI-Virus-SARS-CoV-2 data hub). There are 26 SARS-CoV-2 sequences available from India. We downloaded all Indian SARS-CoV-2 sequences and Wuhan SARS-CoV-2 sequence, which was deposited for the first time after COVID-19 cases started to appear in Wuhan province, China. We downloaded protein sequence of NSP12/RdRp from the database. All these RdRp sequences were first aligned in CLUSTAL Omega to check for similarities or differences. The Clustal Omega algorithm produces a multiple sequence alignment by producing pairwise alignments. Using this algorithm, we identified seven mutations in isolates from Indian SARS-CoV-2 samples compared to Wuhan SARS-CoV-2 RdRp sequence as shown in Fig. 1. We considered the original Wuhan sequence as the wild type for this comparison. These mutations on RdRp are A97V, A185V, I201L, P323L, L329I, A466V and V880I respectively. Out of these, one isolate have triple mutation, two isolates have double mutation, and rest have single mutations (Fig. 1).

P323L mutation causes the alteration in secondary structure of RdRp

Next, we studied the effect of these mutations on secondary structure of RdRp. Our data revealed that mutation at four positions have no effect in secondary structure that includes I201L, L329I, A466V and V880I (Figs. 2C, 2E, 2F and 2G). However, mutation in rest of the three sites causes change in secondary structure (A97V, A185V and P323L) as shown in Figs. 2A, 2B and 2D. At two positions (97 and 185) the alanine amino acid is substituted by valine. The valine side chain is larger than alanine, and substitution of valine at position 97 and 185 impairs packing of the protein as revealed by our secondary structure predications at these two sites (Figs. 2A and 2B, compare i and ii). Our analysis showed that there is addition of two sheet at positions 97 and 98 due to mutation of A97V (Fig. 2A, compare i and ii). Similarly, in A185V mutant, there is a loss of turn at 184 and replacement of helix with sheet structure at 181, 182, 183 and 184 positions (Fig. 2B, compare i and ii). Further, our secondary structure prediction also showed changes in secondary structure when proline is substituted by leucine at 323 position (Figs. 2D and 2H, compare i and ii). The detailed analysis revealed that the mutant RdRp (P323L) have attained considerable changes in secondary structure at the mutation site. There is a loss of turn structures from position 323 and 324 and addition of five sheets at positions 321, 322, 323, 324 and 327 (Figs. 2D and 2H, highlighted in dashed box). Proline possesses a unique property, its side chain cyclizes back on to the backbone amide position, due to which it contributes to the secondary structure formation because of its bulky pyrrolidine ring that places steric constraints on the conformation of the preceding residue in the helix. The substitution of proline to leucine in the mutant RdRp might result in loss of the structural integrity provided by proline. Leucine is a hydrophobic amino acid and generally buried in the folded protein structure. Altogether, the substitution of alanine to valine at position 97, 185 and proline to leucine at position 323 in mutant RdRp is changing the secondary structure of the protein that might have functional consequences.

P323L alters the stability dynamics of tertiary structure of RdRp

To understand the impact of mutations on tertiary structure of RdRp, we performed protein modelling using Dynamut software (Rodrigues, Pires & Ascher, 2018). The Dynamut software provide us the information about the alteration in protein stability and flexibility due to the mutations in the native protein structure. Our data revealed that there is change in vibrational entropy energy (ΔΔS_VibENCoM) between the wild type (Wuhan isolate) and the mutant Indian isolate (Fig. 3A). Vibrational entropy represents an average of the configurational entropies of the protein within single minima of the energy landscape (Goethe, Fita & Rubi, 2015). The negative ΔΔS_VibENCoM of mutant RdRp represents the rigidification of the protein structure and positive ΔΔS_VibENCoM represents gain in flexibility. Here, our data show that the mutation at A185V and P323L lead to rigidification of mutant protein structure (Figs. 3B and 3D). However, the mutation at I201L leads to increase in flexibility (Fig. 3C). Further, we also calculated the free energy differences, ΔΔG, between wild-type and mutant. The free energy differences, ΔΔG, caused by mutation have been correlated with the structural changes, such as changes in packing density, cavity volume and accessible surface area and therefore, it measures effect of mutation on protein stability (Eriksson et al., 1992). In general, a ΔΔG below zero means that the mutation causes destabilisation and above zero represents protein stabilisation. Here, our analysis showed positive ΔΔG for A185V and P323L mutations suggesting that P323L mutation is stabilising protein structure (Fig. 3A); however, we observed negative ΔΔG for I201L mutation indicating its destabilising behaviour (Fig. 3A). The ΔΔG for P323L mutant was 0.908 kcal/mol which is significantly higher than others. Therefore, we also did protein modelling of P323L mutant using two additional structures of RdRp (7BV1, 7BV2) (Yin et al., 2020) and analysed its impact on ΔΔG. Our data show that in both the cases (7BV1 and 7BV2) the P323L mutation is leading to stabilisation of the protein structure with the values of ΔΔG as 0.530 kcal/mol and 0.460 kcal/mol respectively.

We further closely analysed the changes in the intramolecular interactions due to these three mutations in RdRp. Our data showed that it is affecting the interactions of the residues which are present in the close vicinity of alanine, isoleucine and proline. The substitution of wild type residues with mutant residues alters the side chain leading to change of intramolecular bonds in the pocket, where these amino acids resides as shown in Figs. 4A–4C. Therefore, it can be conclusively stated that the three mutations namely, A97V, A185V and P323L, in large, are changing the stability and intramolecular interactions in the protein that might have functional consequences.