Dominant clade-featured SARS-CoV-2 co-occurring mutations reveal plausible epistasis: An in silico based hypothetical model

Abstract

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has evolved into eight fundamental clades with four of these clades (G, GH, GR, and GV) globally prevalent in 2020. To explain plausible epistatic effects of the signature co-occurring mutations of these circulating clades on viral replication and transmission fitness, we proposed a hypothetical model using in silico approach. Molecular docking and dynamics analyses showed the higher infectiousness of a spike mutant through more favorable binding of G₆₁₄ with the elastase-2. RdRp mutation p.P323L significantly increased genome-wide mutations (p < 0.0001), allowing for more flexible RdRp (mutated)-NSP8 interaction that may accelerate replication. Superior RNA stability and structural variation at NSP3:C241T might impact protein, RNA interactions, or both. Another silent 5'-UTR:C241T mutation might affect translational efficiency and viral packaging. These four G-clade-featured co-occurring mutations might increase viral replication. Sentinel GH-clade ORF3a:p.Q57H variants constricted the ion-channel through intertransmembrane-domain interaction of cysteine(C81)-histidine(H57). The GR-clade N:p.RG203-204KR would stabilize RNA interaction by a more flexible and hypo-phosphorylated SR-rich region. GV-clade viruses seemingly gained the evolutionary advantage of the confounding factors; nevertheless, N:p.A220V might modulate RNA binding with no phenotypic effect. Our hypothetical model needs further retrospective and prospective studies to understand detailed molecular events and their relationship to the fitness of SARS-CoV-2.

Keywords: COVID-19; SARS-CoV-2; clades; co-occurring mutations; fitness; infection paradox; virulence.

Figure 1

Schematic diagram of SARS‐CoV‐2 replication in the cell showcasing S, N, ORF3a, RdRp, NSP3, and 5′‐UTR based epistatic interactions. The replication cycle starts with the ACE2 receptor binding of the spike glycoprotein (S) as cornered at the top‐left and finishes with the exocytosis at the top‐right. The viruses which do not carry G‐, GH‐ and/or GR‐featured mutations in the S, N, ORF3a, RdRp, NSP3, and 5′‐UTR are denoted as the wild‐type where mutants contain those. The red and green color icons throughout the diagram, such as proteins, genome, and virion, represent the wild and mutant types, respectively. For a generalized virion, we used the blue color. Although this theme does not show the co‐infection of both types, which might occur on rare occasions, we showed the comparative epistatic effects side‐by‐side fashion during the whole replication cycle, making it easy to grasp. Related figure(s) for each protein are shown in the enclosed box. We used the “question” mark in a pathway and an explanatory box to mean the uncertainty or unknown effects of any mutant proteins/RNA structure. ACE2, angiotensin‐converting enzyme 2; ER, endoplasmic reticulum; ERGIC‐endoplasmic; RdRp, RNA dependent RNA polymerase; NSP14‐proof‐reading enzyme of SARS‐CoV‐2; SARS‐CoV‐2; severe acute respiratory syndrome coronavirus 2; 5′‐UTR, 5′‐untranslated region

Figure 2

Different structural and stability comparisons of the wild and mutant spike protein. Structural superposition of wild and mutant spike proteins (A) and (B); conformation in S1–S2 (C) and S2ʹ sites (d) and (E); representation of vibrational entropy energy change on the mutant type structure (F); and interatomic interaction prediction of both wild (G) and mutant (H) types. For Figure (A)–(E), the gray and yellow colors represent the wild and mutant protein, respectively. (A) The downstream (617–636) of D614G in wild (green) and mutant (red) S protein was focused. Overlapping the wild (D₆₁₄) and mutant (G₆₁₄) S protein showed a conformational change in the 3D structures. (B) However, the conformational change is in the loop region (618–632) of the proteins thus may potentially play a role in interacting with other proteins or enzymes, such as elastase‐2, as we focused on in this study. (C) No change was found in the S1–S2 cleavage site (685–686), depicted in blue color, of the wild and mutant protein. (D) Surface and (E) cartoon (2°) structure of the superimposed wild and mutant proteins where the S2ʹ (pink) is situated in the surface region and does not show any change in accessibility in the residual loop region. (F) The mutant (G₆₁₄) protein showed higher flexibility in the G⁶¹⁴ (sticks) and its surroundings (red). The intramolecular interaction determined the overall stability of the (G) wild and (H) mutant structure where C_β of D⁶¹⁴ (aspartic acid at 614; green stick modeled) had two hydrophobic interactions with the benzene rings. These intramolecular contacts stabilize the S protein of wild‐type. The missing of this bond destabilizes the mutant (G₆₁₄) protein. The mutant protein has glycine at 614, which has less chance of interacting with neighboring aa due to its shorter and nonpolar R‐group. The color code representing the bond type is presented in each (G) and (H)

Figure 3

The molecular docking of wild and mutant with elastase‐2. Both the (upper figure) wild (D₆₁₄) and (lower figure) mutant (G₆₁₄) version of S protein was shown in golden color, whereas the elastase‐2 docked to D₆₁₄ and G₆₁₄ in blue and green color, respectively. The enlarged views of the docked site were shown in separate boxes. (A) The possible docked residues (stick model) on the wild S protein (warm pink) and elastase‐2 (green) are 618(Thr)−619(Glu)−620(Val) and 198(Cys)−199(Phe)−225:227 (Gly, Gly, Cys), respectively. The aspartic acid at 614 is 17.3°A far away from the valine (101), apparently the nearest aa of elastase‐2 to the cleavage site (615–616). (B) The possible interacting residues (stick model) on the mutant S protein (blue) and elastase‐2 (warm pink) are 614(Gly)−618(Thr)−619(Glu)−620(Val) and 101(Val)−103(Leu)−181(Arg)−222:227(Phe, Val, Arg, Gly, Gly, Cys), and 236 (Ala), respectively. In this case, the glycine at 614 is 5.4°Afar from the valine (101), the nearest aa of elastase‐2 to the cleavage site (615–616)

Figure 4

Molecular dynamics of spike‐elastase2 and RdRp‐NSP8 complex (A) Both the wild and mutated spike protein had lower RMSD profile till 60 ns, then it rose and maintained a steady state. Although the spike protein had a higher degree of deviation in the RMSD profile than RdRp, they did not exceed 3.0 Å. The RMSD demonstrated that mutant and wild RdRp protein complex has an initial rise in RMSD profile due to flexibility. Therefore, both RdRp complexes stabilized after 30 ns and maintained a steady peak. The wild‐type RdRp complex had a slightly higher RMSD peak than mutant RdRp, indicating the more flexible nature of the wild‐type. (B) The spike protein complex had a similar SASA profile, did not change its surface volume, and maintained a similar trend during the whole simulation time. The higher deviation of SASA indicates that mutant and wild‐type RdRp had a straight line. Still, the mutant structure had a higher SASA profile, indicating the protein complex had enlarged its surface area. Therefore, the mutation in RdRp protein leads to more surface area expansion than wild types as their average SASA value had a significant difference. (C) Mutated spike exhibits a little more Rg profile than the wild‐type, which correlates with the comparative labile nature of the mutant. The higher level of Rg value defines the loose packaging system and mobile nature of the protein systems. The mutant RdRp had a lower level of fluctuations and maintained its integrity during the whole simulation time. The wild‐type RdRp complexes had higher deviations and more mobility than the mutant complex. (D) Any aberration in hydrogen bond numbers can lead to higher flexibility. Therefore, the mutant and wild spike proteins exhibit the same flexibility in terms of H‐bonding. The mutant RdRp protein had more hydrogen bonding than the wild types, but they did not differentiate too much, and a relatively straight line was observed for the protein. RMSD, root‐mean‐square deviation

Figure 5

The molecular interaction of mutant RdRp with NSP8. The mutant (L₃₂₃) RdRp (pale green) and NSP8 (light blue) are interacting as shown in the center of the lower figure, and an enlarged view of the docked site is presented above within a box. The leucine at 323 interacted with the Asp (112), Cys (114), Val (115), and Pro (116). The wild (P₃₂₃) RdRp has identical docking interactions with NSP8 (Table 1), thus it is not presented as a separate figure

Figure 6

The effect on transmembrane channel pore of ORF3a viroporin due to p.Q57H mutation. (A) The wild (Q₅₇) and mutant (H₅₇) ORF3a proteins are presented in light gray and green colors. The structural superposition displays no overall conformation change; however, the histidine at the position 57 of mutant ORF3a (deep blue) is slightly rotated from glutamine at the exact position of the wild protein (bright red). This change in rotamer state at the residue 57 may influence (B) the overall stability of H₅₇ (upper part) overQ₅₇ (lower part)because of ionic interaction of histidine (green; stick model) of transmembrane domain 1 (TM1) with cysteine at 81 (yellow stick) of TM2. The color code defined different bond types is shown in the inlet

Figure 7

Structural superposition of wild and mutant N protein. The light gray color represents both wild (RG_203–204) and mutant (KR_203–204) N protein. The linker region (LKR: 183–247 aa) of wild (RG_203–204) and mutant (KR_203–204) are in pale yellow and warm pink color, respectively. (A) The aligned structures showed a highly destabilizing (Table 2) conformational change from 231 to 247 aa within LKR. Other regions of the N protein, especially the SR‐rich region (184–204 aa) where the mutations occur, do not change. (B) A more emphasized look into the SR‐rich and mutated sites (RG203–204KR) of wild and mutant N protein represent a slight deviation in the Ser (197) and Thr (198) where only glycine (green) to arginine (blue) substitution at position 204 shows changing at rotamer state. The enlarged view is shown in the bottom part