Emerging mutations in envelope protein of SARS-CoV-2 and their effect on thermodynamic properties.

Structural proteins of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) are potential drug targets due to their role in the virus life cycle. The envelope (E) protein is one of the structural proteins; plays a critical role in virulency. However, the emergence of mutations oftenly leads to drug resistance and may also play a vital role in virus stabilization and evolution. In this study, we aimed to identify mutations in E proteins that affect the protein stability. About 0.3 million complete whole genome sequences were analyzed to screen mutations in E protein. All these mutations were subjected to stability prediction using the DynaMut server. The most common mutations that were detected at the C-terminal domain, Ser68Phe, Pro71Ser, and Leu73Phe, were examined through molecular dynamics (MD) simulations for a 100ns period. The sequence analysis shows the existence of 259 mutations in E protein. Interestingly, 16 of them were detected in the DFLV amino acid (aa) motif (aa72-aa75) that binds the host PALS1 protein. The results of root mean square deviation, fluctuations, radius of gyration, and free energy landscape show that Ser68Phe, Pro71Ser, and Leu73Phe are exhibiting a more stabilizing effect. However, a more comprehensive experimental study may be required to see the effect on virus pathogenicity. Potential antiviral drugs, and vaccines may be developed used after screening the genomic variations for better management of SARS-CoV-2 infections.

Introduction

Recently the newly emerged n class="Disease">coronavirus disease-19 (COVID-19) remains a major public health issue since 2019. The causative agent severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) [1] is passing through numerous evolutionary stages [2,3] and seems more fatal. The virus is rapidly spreading through respiratory droplets, affecting the important organs of the host body [2,3].

Since the coronavirus 2 (n class="Species">CoV-2) belongs to the formerly known family of coronaviruses, it shares close genomic similarity with SARS-CoV. The CoV-2 contains the ssRNA genome [6]; it encodes four structural proteins such as a spike (S), an envelope (E), membrane (M), and nucleocapsid (N) [7], resulting in the coding of 16 non-structural proteins (NSPs). These structural proteins are responsible for viral replication and virion-receptor attachments and, thus, are involved in pathogenicity, viral spread, and the introduction of the virus to the cells of the host body.

The E protein is con class="Gene">mposed of 76 amino acids [8] weighing 8–12 kDa [9]. This protein has an N-terminal domain (NTD), hydrophobic domain, and C-terminal domain (CTD) [10], arranged as NTD first 1–9 amino acid (aa), the hydrophobic area, extending from 10 to 37 aa and CTD (38–75 aa) [11]. Ionic pores are formed across the membrane due to the arrangement of the hydrophobic tail region. Structurally, the E protein consists of seven alpha helices and eight loop regions. It forms a pentameric configuration (five molecules) with 35 alpha-helical regions and 40 looped ring regions, which are formed by the hydrophobic tail (Fig. 1 B). This pentameric configuration of the hydrophobic tail can be affected by the interactions within CTD [12]. These pores serve as ion channels that allow the movement of the virus across the membrane, enhancing its pathogenicity [13]. E mutations reduce viral pathogenicity and also stop the channel activity [14] representing an essential drug target and vaccine candidate [15].

Fig. 1

Domain organization of E proteins and the location and frequency of most common mutations. TM: transmembrane [15], NTD: N-terminal domain, CTD: C-terminal domain. (B). Pentameric representation of E proteins showing channel (C). Full-length E protein (I-TASSER, E-QHD43418) with most common mutations in the loop region.

Domain orgaene">nizatioene">n of n class="Gene">E proteins and the location and frequency of most common mutations. TM: transmembrane [15], NTD: N-terminal domain, CTD: C-terminal domain. (B). Pentameric representation of E proteins showing channel (C). Full-length E protein (I-TASSER, E-QHD43418) with most common mutations in the loop region.

Some novel n class="Gene">mutations seen in the E protein sequence do not present in the already existing coronaviruses. In the CoV-2 E protein sequence, arginine is mutated by isoleucine, threonine, and lysine (R69I, R69T, and R69K) at the 69th position. Moreover, in the amino acid sequence of CoV-2, serine and phenylalanine are present at 55th and 56th positions, instead of threonine and valine, respectively [16].

The SARS-CoV-2 has undergoene">ne nun class="Gene">merous mutations in all the important targets [[17], [18], [19]]; therefore, drugs created to fight CoV-2 might not very potent. In order to evaluate the drug targets for designing novel antiviral drugs against CoV-2, it is vital to screen the frequency of mutations and their effect on the thermodynamic properties of the important target proteins. This is particularly important for the effective treatment of emerging microbes, including mutants of SARS-CoV-2.

Although structural proteins of CoV-2 are in class="Gene">mportant for investigating mutations, here in this study, we investigated the variations existing only in E proteins owing to their potential drug target. Very limited information and small-scale genomic data has been screened for mutations in the E protein. This is the first comprehensive study in which we screened 0.295 million complete genomes of SARS-CoV-2 for identifying variants in the E protein. Among all the genomes, 259 mutations were detected, exhibiting various degrees of thermodynamic properties.

Analyzing the frequency of mutatioene">ns aene">nd their effect oene">n the thern class="Gene">modynamics properties of E proteins may allow the researcher to design novel inhibitors and predict the level of pathogenicity and transmission. In the current study, the most common mutations that were detected at the C-terminal (Ser68Phe, Pro71Ser, Leu73Phe) were exhibiting more stabilizing effect on E proteins structure. First, the presence of a large number of mutations in the E protein of CoV-2 may lead to the conformational changes and evolutions, resulting in therapeutic failure. Secondly, the virus may be more fatal in near future, and therefore, antivirals against geography-specific CoV-2 strains might be more effective.

Methodology

Genomic sequence retrieval

The complete genon class="Gene">me was retrieved from the Global Science Initiative on Sharing All Influenza Data (GISAID) (December 2019–2020) (https://www.gisaid.org/) [20]. The GISAID shared the virus data to publish results and the metadata relevant to public health scientists. This server provides all kinds of CoV-2 genomic data, including even those that have not been unpublished. We screened 0.295 genomes of CoV-2 reported worldwide for variants analysis in E proteins. The sequences were aligned with the reference CoV-2 genome (Accession NC_045512) using the CoVsurver application (https://www.gisaid.org/epiflu-applications/covsurver-mutations-app/). The identified mutations in structural proteins of CoV-2 were separated and arranged in the form of excel sheets. The statistical analysis was performed to screen the most common variants.

Structural information

Scientists are all well aware that sharing the genomic aene">nd proteon class="Gene">mic data of SARS-CoV-2 is important for the better management of infectious diseases in order to devise countermeasures. The E protein structural data was retrieved from the protein data bank (PDB) [21] (PDB IDs: E = 7k3g). Some residues at the NTD (1–7aa) and CTD (39–75) are missing in the PDB structure of the COV-2 E protein. To observe the effect of mutations on the residues missing in the CTD and ND terminals, the full-length E protein structure was downloaded (E-QHD43418) from I-TASSER [22,23]. The chain ID in the I-TASSER structure is missing, which was added in PYMOL using the code “sele, chain'” and “alter (sele), chain = ‘A'.” The I-TASSER has already modeled the full-length proteins of CoV-2 using the NCBI reference data (NC_045512) (GenBank MN908947).

Mutation effect on E proteins’ dynamic stability

All the observed mutatioene">ns were recorded aene">nd their effect oene">n the n class="Gene">E protein structure was computed using the DynaMut [24]. The server implemented the mutation effect and also the normal mode methods that can be used to analyze the variants that affect protein stability and flexibility. This impact is measured through graph-based signatures and also as a normal mode. This method outperforming (P < 0.001) with results are also displayed in good resolution.

Molecular dynamics (MD) simulation

The MD sin class="Gene">mulation was performed on the Desmond module (Schrodinger), as described in the previous study [23,24]. Briefly, the TIP3P model and Gromos9643a1 forefield was applied. The system was neutralized with counterions (NaCl). The cubic simple point charge (SPC) water box was applied. Two-step (NVT and NPT) energy minimization (50000 ps) till the minimization completion was continued. The ambient pressure was set at 1.013 bar and temperature 310 K for 100 ns? The thermodynamic stability of the WT (wild type) and mutants (MT) E proteins was analyzed using root means square fluctuations (RMSF) and root means square deviation (RMSD). All the simulation were repeated three time for better results. Approximately 1000 frame per simulation was run.

The Gibbs free energy

The Gibbs free energy (G) [27] of MTs was plotted against wild type n class="Gene">E protein. The G is minimized to equilibrium state of system at constant temperature and pressure which is a thermodynamic potential. Principal component analysis (PCA) is performed to recognize low modes in proteins [28,29]. PCA simplifies the complicated motion in trajectory [[30], [31], [32]]. A set z1, z2 …, zp known as principal components (PCs) were generated during PCA. Energies of sets of proteins conformations is called Free Energy Landscape (FEL) [33,34]. The first two components (PC1 and PC2) give the trajectories on initial two principal components of motion. G values shows the stability level of proteins [[35], [36], [37]].

Results and discussion

This is the first comprehensive study in which 0.295 n class="Gene">million complete genomic sequences of CoV-2, which were reported worldwide in the GISAID server (from December 2019 to December 2020), were analyzed to identify variants in E proteins. A large number of non-synonymous mutations (259) were detected in the genome sequences from 48 countries (Supplementary Table S1), among which the largest number were present in isolates from England (S1). This wide range of variations may project the variation level of CoV-2 strains worldwide. Previous studies [25,26] have screened 3617 and 81,818 COV-2 genomes for E protein variants, respectively. These recent studies reported 115 and 15 non-synonymous mutations, respectively, which were mainly present in the CTD.

Mutations in envelope (E) protein

Mutations in the NTD

Being the smallest structural protein (75aa) of n class="Species">CoV-2, all of the residues’ positions of the E protein harbored non-synonymous mutations (S1). A total of 31 non-synonymous mutations were detected in the NTD (Supplementary File S1), among which V5F (n = 39), E8D (n = 35), V5I (n = 17), and Y2H (n = 14) were common. A majority of these mutations were detected in the isolates from the UK. The NTD helps the tail region target the Golgi complex using some of its associated elements as mutations in this region may affect its efficiency.

Mutations in the CTD

The most frequent n class="Gene">mutations were detected in the CTD of the E protein (68–73 aa) (Table 1 ). Some of them were T91, S55F, V62F, S68F, R691, P71L, P71S, and L73F (Fig. 1). Other common mutations have been listed in the table along with their residue types and positions (Table 1 and S1). Mutations in the CTD, namely, S55F (128), V62F (129), and R69I (159) may affect the virus pathogenesis, altering the binding of the E protein to a tight junction. The motif “DLLV” (72–75 aa) of the CTD showed some mutations that may affect the PALS1 (Protein Associated with Caenorhabditis elegans Lin-7 protein 1) at the Golgi complex, as well as CoV-2 infectivity [38,40]. The PALS1 is a member of the post-synaptic density protein-95/Discs Large/Zonula occludens-1 (PDZ) domain-containing proteins group that is involved in diverse cellular functions and also functions as scaffolds for signaling protein [28,29].

Table 1

Frequency of some common mutations in the CoV-2 E protein.

Accession	Wild type AA	Position	Mutated AA	Frequency	aMutation
EPI_ISL_476911	T	9	I	168	T9I
EPI_ISL_424214	S	55	F	128	S55F
EPI_ISL_538676	V	62	F	129	V62F
EPI_ISL_448073	S	68	F	b419	bS68F
EPI_ISL_452908	R	69	I	159	R69I
EPI_ISL_577907	P	71	L	158	P71L
EPI_ISL_660339	P	71	S	b264	bP71S
EPI_ISL_478788	L	73	F	b218	bL73F

Full list is available in the S1 supplementary file.

; mutations subjected to MD simulations.

Frequency of some con class="Gene">mmon mutations in the CoV-2 E protein.

; mutatioene">ns subjected to n class="Gene">MD simulations.

Mutations in the transmembrane (TM) domain

Transmembrane variaene">nts such T9I (n = 168), n class="Mutation">F20L (n = 90), L21F (n = 84), V24 M (n = 76), and T30I (n = 72) may affect the homo pentameric configuration of the E protein [43]. The E protein may also be an effective drug target as it contributes equally to the pathogenicity and cytotoxicity of the virus. It produces viroporins that are hydrophobic in nature [44]. Proline residues facilitate the targeting of the cis-Golgi complex by the hydrophobic tail present in the cytoplasm. The release of these virion particles is facilitated by the ionic gradient present in the endoplasmic reticulum and Golgi compartment through the E protein [8]. Studies on E mutants for highlighting the structural changes of the ion-channel activity behind mutations, might be very helpful for better management of COVID-19.

The CTD of E proteins harbored son class="Gene">me common mutations whose stability was predicted. The DynaMut prediction outcome of L73F (ΔΔG: −0.417 kcal/mol), P71S (ΔΔG: −0.255 kcal/mol) exert a destabilizing effect. However, T9I (ΔΔG: 0.190 kcal/mol), P71L (ΔΔG: 0.012 kcal/mol), and S68F (ΔΔG: 0.362 kcal/mol) shows a stabilizing effect (Supplementary File S1). The most common mutations were detected at the C-terminal (Ser68Phe, Pro71Ser, and Leu73Phe) were also assessed through MD simulations, exerting a stabilizing effect on the E protein of CoV-2.

Thermodynamic properties

We analyzed the thermodyene">nan class="Gene">mic properties of most common variants (Ser68Phe, Pro71Ser, and Leu73Phe) present in the CTD in relation to the E structure stability in comparison with the wild type (WT).

RMSD and RMSF of WT and MTs

The RMSD graphs of then WT aene">nd n class="Gene">MTs E proteins are shown in Fig. 2 . The MTs S68F and P71S seem more stable from 55ns to 100 ns, exhibiting 9 Å and 10.5 Å RMSDs, respectively. The RMSD of the WT still exhibiting fluctuations at 100 ns? Similar to the other MTs, L73F is also exhibiting stable deviations through the simulation period, with 13.7 Å to 14.5 Å RMSD. The WT exhibited 5.3 Å RMSD at the 80ns–95ns period and still seems to rise from 95ns to 100 ns?

Fig. 2

RMSD comparison of the WT and MTs E proteins of SARS-CoV-2.

RMSD con class="Gene">mparison of the WT and MTs E proteins of SARS-CoV-2.

The comparisoene">n of the WT aene">nd n class="Gene">MTs in terms of flexibility exhibited significant variations at some positions (Fig. 3, Fig. 4 ). The RMSF at the 53rd to 75th aa position exhibited a significant difference between the WT and MTs. The WTs at this position exhibited 6.89 Å (55 aa), 9.88 Å (71 aa), and 9.5 Å (75 aa); S68F demonstrated very low RMSF, below 4 Å, for the majority of the residues and 5.1 Å at 75 aa, which was significantly low as compared to the WT. Similarly, P71S also exhibited very low RMSF at CTD (5.7 Å) when compared to the WT (Fig. 3). Similar to S68F, the MT L73F also demonstrated low-level fluctuation when compared with the WT. MD simulations explore the insight dynamic changes at the molecular level [[45], [46], [47]], which might be difficult to accomplish through experimental work. Several studies have reported that any change in protein functions might be associated with RMSF [[48], [49], [50]]. Flexibility is one of the key thermodynamic characteristics that maintain the optimal functions of proteins [51]. A large change in this property may alter the function of the biomolecules.

Fig. 3

The RMSF of the WT and MTs in the E proteins of SARS-COV-2. CA: Carbon alpha. The RMSF of WT (orange) and MTs. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 4

The RMSDs of WT and MTs at different simulation periods. The WT behaves very differently at the CTD.

The RMSF of the WT aene">nd n class="Gene">MTs in the E proteins of SARS-COV-2. CA: Carbon alpha. The RMSF of WT (orange) and MTs. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

The RMSDs of WT aene">nd n class="Gene">MTs at different simulation periods. The WT behaves very differently at the CTD.

Protein secondary structure elements (SSE) such as alpha helices aene">nd beta straene">nds are n class="Gene">monitored throughout the simulation. The difference in the SSE of the WT and MTs has been highlighted in Fig. 5 at 20ns–70 ns? The mutation induced some changes in the SSE within a 100ns simulation period. Upon examining Fig. 5, some differences in the SSE of the WT (Fig. 5A) and MTs (Fig. 5B, C, D) can be seen.

Fig. 5

Comparison of the SSE of the WT and MTs in MD simulation. The plot reports SSE distribution by residue index throughout the protein structure and summarizes the SSE composition for each trajectory frame over the course of the simulation. The plot monitors each residue and its SSE assignment over a certain period of time. (A) The WT residue index and SSE along Y-axis shows differences at the 60–65ns MD simulations period; (B) S68F shows a little variation at 40–55ns; (C) P71S exhibits more variations at 45–68ns; (D) L73F also demonstrates variations when the WT is compared with P71S and S68F.

Comparisoene">n of the SSE of the WT aene">nd n class="Gene">MTs in MD simulation. The plot reports SSE distribution by residue index throughout the protein structure and summarizes the SSE composition for each trajectory frame over the course of the simulation. The plot monitors each residue and its SSE assignment over a certain period of time. (A) The WT residue index and SSE along Y-axis shows differences at the 60–65ns MD simulations period; (B) S68F shows a little variation at 40–55ns; (C) P71S exhibits more variations at 45–68ns; (D) L73F also demonstrates variations when the WT is compared with P71S and S68F.

Rg of WT and MTs

The degree of protein-folding stability could be measured through Rg. Fluctuatioene">ns in Rg for a period indicate unstable folding, while a straight value reveals stable folding [[52], [53], [54], [55]]. A protein with n class="Gene">misfolding shows variations in Rg over time (Fig. 6 ). The WT and MTs exhibited variations in folding. The MT S68F exhibited more stable folding than the WT, whereas P71S and L73F also presented stable folding for the 18.9–50ns period (Fig. 6). The S68F MT shows the lowest fluctuations from 21ns to 100 ns (1.3 nm). Rg is a mass-weight root mean square distance between atoms and their common center of mass. It is a vital parameter for defining dynamic stability, offering an insight mechanism of dimension and compactness of biomolecules and total protein systems. The Rg plotted for three MTs of the CoV-2 E proteins shows notable stability than that of the WT. The mutant structures show a slightly low average Rg value than the WT.

Fig. 6

Comparison of Rg of the WT and MTs E proteins. The WT (blue) exhibited a little difference in folding throughout the 100ns simulation period as compared to the MTs, showing that the folding in both types is stable. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Comparisoene">n of Rg of the WT aene">nd n class="Gene">MTs E proteins. The WT (blue) exhibited a little difference in folding throughout the 100ns simulation period as compared to the MTs, showing that the folding in both types is stable. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Hydrogen bonding of the WT and the MTs

Most of the interactioene">ns an class="Gene">mong the atoms are mainly hydrogen bonding (HB), which assists in protein folding, stability, and is also involved in recognition. The α helix and β sheet stabilize by the HB potential between the amide nitrogen of the protein and the carbonyl oxygen of the main chain [[56], [57], [58]]. The primary forces in protein–drug interactions are HB, van der Waals, and electrostatic forces [51]. The average number of HB among the WT and MTs CoV-2 E protein atoms is shown in Fig. 7 . MT S68F exhibited fewer HB than P71S and P73F. However, both the MTs (P71S and P73F) demonstrated more HB than the WT (Fig. 7), signifying the stabilizing effect of the mutation. P71S exhibited maximum HB (295) during the first 25ns simulations. However, the number of HB in the MD simulation in the last 25ns is almost the same for the WT and MTs.

Fig. 7

Hydrogen bonding of the WT and MTs E proteins.

Hydrogen boene">nding of the WT aene">nd n class="Gene">MTs E proteins.

The effect of mutatioene">ns oene">n n class="Gene">E proteins’ stability and flexibility has been shown in Fig. 8 . All of the MTs exhibited an increase in flexibility, which may have a positive effect on function. The molecular activities are linked with the flexible regions. According to a previous investigation [59], protein requires flexibility to engage in good catalytic activity. Adding water improves flexibility, and understanding the diverse role of protein flexibility can help develop biotechnological solutions including vaccines.

Fig. 8

WT and MTs stability and flexibility.

WT and MTs stability aene">nd flexibility.

MTs showed n class="Gene">more interactions than WT. The S68F exhibited a different geometry due to the aromatic nature of phenylalanine substitutions. Moreover, the number of interactions in this MT seems greater than in the WT (Fig. 9 ). In MT P71S, the proline has been substituted into serine, present in the active site of many enzymes. Similar to the S68F, the L73F exhibited more interactions (Fig. 9) with surrounding residues than the WT due to its substitution into aromatic amino acid, forming more hydrogen interactions with its surrounding residues that might be instrumental in E stability.

Fig. 9

Interactions of WT and MTs residues with surrounding aa. The type of interactions of WT and MTs residues is colored-coded. S68F, P71S, and L73F MTs have been compared with their WT on the left.

Interactions of WT and MTs residues with surrounding aa. The type of interactioene">ns of WT aene">nd n class="Gene">MTs residues is colored-coded. S68F, P71S, and L73F MTs have been compared with their WT on the left.

The FEL of WT and MTs has been showene">n in Fig. 10 . The lowest aene">nd stable energy state is represented by blue color. WT n class="Gene">E protein showed more stable state that two MTs (P71S, L73F) whereas the MT S68F seems more stable than WT. The blue areas show stability while other indicate transitions in the protein conformation to attain a more favorable state. The green areas in P71S, L73F are more prevalent that WT which shows stability next to blue regions. These result demonstrates that E protein MTs might be useful to cause viral pathogenecity [[60], [61], [62]]. Calculating G might be important to observe the overall stability upon mutations.

Fig. 10

Gibbs free energy landscape. The scale shows the free energy values. Blue and green regions are more stable than red and yellow. MT S68F shows more stability than WT based on the free energy landscape. WT exhibited more stability than P71S and L73F. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Gibbs free energy landscape. The scale shows the free energy values. Blue and green regions are more stable thaene">n red aene">nd yellow. n class="Gene">MT S68F shows more stability than WT based on the free energy landscape. WT exhibited more stability than P71S and L73F. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Conclusion

In this comprehensive study, 259 noene">n-syene">noene">nyn class="Gene">mous mutations were detected in the E protein with different frequencies. Mutations were present in all three domains—the CTD, the NTD, and the hydrophobic domain; however, the highest frequency was detected in the first one. All of the MTs demonstrated various degrees of flexibilities and stabilities, where the majority depicted a loss in flexibility and stability. Moreover, three MTs (S68F, P71S, and L73F) were analyzed through MD simulations and exhibited a stabilizing effect on the E protein structure. In order to evaluate the drug targets for designing a novel antiviral drug against CoV-2, it is vital to dig out the frequency of mutations and their effect on the thermodynamic properties of the selected target. This is particularly imperative for the effective treatment of emerging microbes including SARS-CoV-2 mutants. The current study will contribute significantly to the knowledge on virus stabilization and the evolutionary aspects and pathogenicity of CoV-2 infections, which might be useful for better management of COVID-19 and the development of a vaccine in the near future.

Declaration of competing interest

All the authors have no competing interests.