Unraveling the Molecular Mechanism of Recognition of Human Interferon-Stimulated Gene Product 15 by Coronavirus Papain-Like Proteases: A Multiscale Simulation Study.

The papain-like protease (PLpro) of the coronavirus (CoV) family plays an essential role in processing the viral polyprotein and immune evasion. Additional proteolytic activities of PLpro include deubiquitination and deISGylation, which can reverse the post-translational modification of cellular proteins conjugated with ubiquitin or (Ub) or Ub-like interferon-stimulated gene product 15 (ISG15). These activities regulate innate immune responses against viral infection. Thus, PLpro is a potential antiviral target. Here, we have described the structural and energetic basis of recognition of PLpro by the human ISG15 protein (hISG15) using atomistic molecular dynamics simulation across the CoV family, i.e., MERS-CoV (MCoV), SARS-CoV (SCoV), and SARS-CoV-2 (SCoV2). The cumulative simulation length for all trajectories was 32.0 μs. In the absence of the complete crystal structure of complexes, protein-protein docking was used. A mutation (R167E) was introduced across all three PLpro to study the effect of mutation on the protein-protein binding. Our study reveals that the apo-ISG15 protein remains closed while it adopts an open conformation when bound to PLpro, although the degree of openness varies across the CoV family. The binding free energy analysis suggests that hISG15 binds more strongly with SCoV2-PLpro compared to SCoV or MCoV. The intermolecular electrostatic interaction drives the hISG15-PLpro complexation. Our study showed that SCoV or MCoV-PLpro binds more strongly with the C-domain of hISG15, while SCoV2-PLpro binds more favorably the N-domain of hISG15. Overall, our study explains the molecular basis of differential deISGylating activities of PLpro among the CoV family and the specificity of SCoV2-PLpro toward hISG15.

Introduction

The current worldwide outbreak by the severe acute respiratory syndrome-related coronavirus 2 (SARS-CoV-2 or SCoV2) caused the spread of the disease (COVID-19) on an unprecedented global scale.[1] It has emerged as a highly contagious novel coronavirus (CoV).[2,3] According to the World Health Organization (WHO), more than 248 million people have been infected worldwide, including more than 4.9 million deaths, as of November 3, 2021. Consequently, WHO declared the situation a global health emergency. The novel CoV is considered far more contagious and catastrophic than other flu viruses, including SARS-CoV (SCoV) and MERS-CoV or MCoV (Middle East respiratory syndrome coronavirus).[4] SCoV and MCoV caused the outbreak in 2003 and 2012, respectively. The scientific community worldwide is making significant efforts to explore diverse mechanisms targeting viral replication via repurposing approved drugs using experimental and computational approaches.[5,6] However, to date, no effective therapeutic agent has been approved against COVID-19 except a few tested drugs like chloroquine derivatives, remdesivir, azithromycin, favipiravir, ivermectin, etc.[7−10]

SCoV2 is a positive-stranded RNA Betacoronavirus.[11] Usually, these viruses harvest polypeptides that promote proteolytic cleavage in their life cycle. Two crucial proteases, namely the main protease (Mpro) and papain-like protease (PLpro), make the functional replicase complex, leading to viral spread.[12] Meanwhile, various efforts have been made to explore these proteases in combating this noxious COVID-19.[13−15] Great attention is given to Mpro because of its post-translational modification activity of replicase polyproteins.[14] It comprises three domains and a catalytic dyad (Cys145 and His41). On the other hand, PLpro includes a catalytic triad (C111, H272, D286 in SCoV2). The PLpro of SCoV2 and SCoV share ∼83% similarity, while only ∼31% similarity is observed between PLpro of SCoV2 and MCoV (see Figure A). Apart from cleaving polyprotein, PLpro has additional deubiquitination and deISGylating activities that promote virus replication. Deubiquitinating enzymes (DUBs) can remove the post-translational modification Ub from target proteins.[16] Similarly, deISGylating enzymes reverse the post-translational modification of the Ub-like (UbI) protein interferon-stimulated gene product 15 (ISG15) from cellular proteins.[17] The SCoV2/PLpro subdomains comprise UBI, thumb, finger, and palm (see Figure B). The flexible BL2 loop recognizes the substrate. It undergoes large conformational changes from close (unliganded) to open (liganded).[18] During viral infection, the first line of defense against viral pathogens is the innate response, involving the elevated expression of IFN stimulated genes (ISG15) via IFN-cytokine signaling pathways. PLpro cleaves ubiquitin and ISG15 from the host cell proteins that aid CoVs to circumvent the host innate immune responses.[19] The dual functionality of PLpro makes it an attractive antiviral drug target.[20−22]

Figure 1

(A) Multiple sequence alignment (MSA) of PLpro’s of MERS-CoV, SARS-CoV, and SARS-CoV2 using Clustal Omega. Meaning of different symbols used in the MSA are: “*” indicates perfect alignment, “:” indicates a site belonging to a group exhibiting strong similarity, and “.” indicates a site belonging to a group exhibiting weak similarity; (B) the cartoon representation of the PLpro of SCoV2, showing its different regions; thumb (green), finger (blue), palm (pink), UBI domain (brown), and BL2 loop (red);(C) the cartoon representation of the full-length hISG15, showing its N- (purple) and C-domain (cyan) and the connecting hinge region (gold);(D) the complex structure of the SCoV2-PLpro’s (orange) with the full-length hISG15 (sea green).

ISG15 is an interferon-regulated gene that regulates various cellular pathways involved in the host immune responses against different viral infections.[23−25] ISG15 contains two ubiquitin-like domains connected by a hinge or linker, making a linear dimeric protein[26] (see Figure C). Relative to ubiquitin, ISG15 has low cross-species conservation, ranging from 58% to as low as 35% in mammals.[26] Overall, CoVs-PLpro has been reported with robust deubiquitinating activities.[17,27] However, few studies have unveiled its importance in deISGylating activities.[28,29] The crystal structures of PLpro with the C-domain of ISG15 have been reported for all major coronavirues.[30−32] However, it is important to consider the full-length hISG15 since both domains of hISG15 interact with PLpro.[33]

Herein, we study the structural dynamics of the full-length hISG15 and its interaction with the PLpro of SCoV, MCoV, and SCoV2 via atomistic molecular dynamics (MD) simulations in conjunction with the molecular mechanics/Poisson-Boltzmann surface area (MM/PBSA) scheme. First, the complex structure of full-length hISG15 with CoVs-PLpro was prepared and analyzed at the molecular level, followed by MD simulations. In addition, to elucidate their molecular recognition and specificity among CoVs, we performed free energy calculations. Overall, the current study provides insights into the molecular basis of CoVs-PLpro deISGylating activities, especially for the novel SARS-CoV-2 specificity toward hISG15 that may be exploited in targeting SARS-CoV2-PLpro to combat COVID-19. To our knowledge, this is the first study where the structural dynamics of the full-length hISG15 protein against all the major CoVs were explored in detail.

Materials and Method

System Preparation

The complex structure of PLpro with the C-terminal domain of hISG15 for MCoV and SCoV was obtained from Protein Data Bank (PDB) with an accession code of 5W8U[30] and 5TL6,[32] respectively. However, the crystal structure of the full-length Apo-hISG15 is already available (PDB ID: 1Z2M[34]), and the same was downloaded. Next, we superimposed the coordinates of the full-length Apo-hISG15 with the available complexes of MCoV and SCoV to obtain the final complex structures. In the case of SCoV2, the structure of PLpro was obtained from the PDB database (PDB ID 7JRN[31]), which was further superimposed with the constructed SCoV-PLpro-hISG15 complex to model the SCoV2-PLpro-hISG15 complex. (Figure D). We also estimated the structural differences between our construct and the available partial crystal structure by superimposition, which yields no significant differences. The root mean square deviation (RMSD) difference is below 0.5 Å, suggesting a good initial complex structure to run the MD simulation.

Further, we considered one charge-flip mutation R167E that was reported to be 20 times less efficient in hydrolyzing ISG15 by SCoV-PLpro.[32] Hence, we studied the effect of the same mutation in SCoV2, which closely resembles SCoV. Also, the impact of the mutation was investigated for the diverse family of CoV, i.e., MCoV.

In the current study, six complex systems of PLpro-hISG15 were prepared for MD simulations: three wildtypes (WTs) (SCoV, SCoV2, and MCoV) and three mutants (SCoV(MT), SCoV2(MT), and MCoV(MT)). Besides, in order to provide structural dynamics of hISG15 upon complexation with PLpro, we performed MD simulations of Apo-hISG15.

Simulation Protocol

All simulations were conducted using the pmemd.cuda module of AMBER18.[35] Missing hydrogens in the crystal structures were added using the LEaP module of Ambertools19 followed by neutralizing all systems by adding a proper amount of counterions. Prior to the system preparation, protonation states of the charged residues were determined using the Propka 3.1 webserver.[36] All systems were then solvated into periodic octahedron TIP3P water box,[37] keeping a buffering distance of 10 Å from all sides. The Amber ff14SB force field[38] was used for both Apo and complex forms for the simulations. A 2.0 fs time step was fixed for all these seven simulation systems. The SHAKE algorithm[39] was used to constraint the bond, including hydrogen atoms. The particle-mesh Ewald summation (PME)[40] was used to compute long-range interaction with a 10 Å nonbonded cut-off. At first, systems were subjected to energy minimization using 500 steps of steepest descent followed by 500 steps of the conjugant gradient algorithm. Amino acids were kept fixed in this step using the force constant of 2.0 kcal mol–1 Å–2. Secondly, in the absence of restraint, the entire systems were minimized on the solute using the steepest descent algorithm followed by the conjugant gradient algorithm. The systems were then gradually heated from 0 to 300 K in the NVT ensemble, where protein atoms were fixed using a force constant of 2.0 kcal mol–1 Å–2. The temperature was maintained using a Langevin thermostat.[41] The equilibration was conducted for 1 ns in the NPT ensemble without any restraint on the systems. Finally, both wild type and mutant complexes underwent a 2 × 2 μs production run. For the Apo-hISG15 system, the production simulation was conducted for 2 × 4 μs. Overall, 200,000 and 400,000 snapshots were generated for complexes and Apo, respectively.

All analyses, such as RMSD, solvent accessible surface area (SASA), the radius of gyration (Rg), and other analyses were performed using the Cpptraj module of AmberTools19.[42] For hydrogen-bond calculations, a distance cut-off ≤3.5 Å and an angle cut-off of ≥120° were used.[43]

Principal Component Analysis

The principal component analysis (PCA) scheme has been applied to identify the lowest energy states in the collected structural and energetics MD trajectory data.[44] It is a dimensional reduction method and we were able to extract the dominant modes from the whole trajectory. In this technique, the covariance matrix obtained from the fluctuation of amino acid residues is diagonalized, and their motion is represented in terms of their eigenvectors and eigenvalues. The eigenvectors represent the direction of motions, while their respective eigenvalues give the amplitude of motions. The original data projected on the eigenvectors yield the principal components (PCs). The covariance matrix was calculated using the eq :where X and X are the Cartesian atomic coordinates of Cα atoms i and j, while and denote the mean of ith or jth atoms over the ensemble. Also, to investigate the dynamics in greater detail, the hinge region was subjected to dihedral PCA (dPCA), which uses the dihedral angles of each peptide bond of a selected region of interest.

Further, we also generated a free energy surface using the first to PC components for both the PC analysis, a common technique to identify the potential barrier between the different conformational states. This was calculated using the formula ΔG = kBTln(ρ), where kB is the Boltzmann constant, T is the temperature, and ρx is the probability density of the geometric coordinate x.

Binding Free Energy Analysis

The effective interaction between PLpro and hISG15 in terms of free energy values was calculated using the Molecular Mechanics Poisson-Boltzmann Surface Area (MM/PBSA) scheme.[45−49] The binding free energy (ΔGbind) was estimated using the MMPBSA.py script of the AMBER package. This method is less accurate compared to other computationally expensive methods like free energy perturbation and thermodynamic integration. However, this scheme offers a good compromise between speed and accuracy and yields good qualitative agreements to the experimental observations.[50] The MM-PBSA protocol was already discussed in detail in our previous studies.[51−54] Here we used the same protocol for calculating the binding free energy. The following equations estimate the binding free energy in the MM/PBSA scheme:[55]

where ΔEinternal, ΔGsolv, and TΔS represent total internal energy, desolvation free energy, and conformational entropy, respectively. The total internal energy (ΔEinternal) still further consists of ΔEcov (bond, dihedral, and angle), ΔEelec (electrostatic), and ΔE (van der Waals) and is represented by the following equation,

while the desolvation free energy composed of polar solvation (ΔG) and the nonpolar solvation energy (ΔG) is represented by,

The polar part of the solvation energy ΔGpol was calculated using the linearized Poisson-Boltzmann model with the dielectric constant of 1.0 and 80.0 for solutes and solvent, respectively, while the nonpolar part was computed using the surface area. A total of 10,000 conformations was used to estimate the binding free energy from the last 1 μs timescale in the equal interval. The entropic contribution was not considered due to the high computational cost.

Further, at the residual level, the decomposition of the total binding free energy, which comprises van der Waals, electrostatic, polar, and nonpolar terms, was conducted by the Molecular Mechanics Generalized-Born Surface Area method.

Residual Network Analysis of Protein

The protein contact network or residual interaction network can represent the protein structure in the form of networks. In the present study, the network analysis of protein structure (NAPS) server was used to construct protein networks.[56] The last snapshots from the production simulations were used as input. The network can be drawn based on Cα, Cβ, or centroids. Here, nodes in a network denote the Cα atom pairs, and edges between Cα-Cα residue pairs are drawn based on the threshold distance of ∼7 Å. Various network-based hydrophobic, hydrophilic, or charged contacts were analyzed to identify the essential residual interactions in the network. The long-range interaction in a protein network can be found using the shortest-path parameters. The Floyd–Warshall[57] algorithm was used to calculate the shortest path, which defines the minimum number of nodes required to reach from one node to the other.[58,59]

Results and Discussion

Apart from the proteolytic activity, PLpro has deubiquitinating and deISGylating activities that make it a potential antiviral target. Hence to understand and explore the underlying interaction mechanism between PLpro and hISG15, we have provided an atomistic insight into the structural and binding free energy analysis using the large-scale MD and /MM-PBSA scheme. Moreover, in-depth analysis of their molecular recognition will help design inhibitors that could inhibit the interaction between SCoV2-PLpro and hISG15.

Structural Dynamics of Systems

The convergence of simulations can be noted from the time evolution of the RMSDs of protein’s backbone atoms, as shown in Supporting Information Figures S1–S3. Further, to confer the conformational dynamics states of hISG15 in its Apo and complex form, both in the WT and the mutant complexes of CoVs-PLpro, the combined probability distributions of the respective RMSDs are shown in Figure CA. The distinctive peaks from the plots correspond to the most probable conformational states of a system. Figure A reveals that the charge-flip mutation (R167E) systems exhibited many distinct conformations at varying RMSD values than their respective WT, like for the SCoV2(MT) (at ∼5 and ∼9 Å) and MCoV(MT) (at ∼7.3 and ∼8.2 Å) complexes. The SCoV complex showed a quite high stable conformation at ∼5 Å, and its mutant showed nearly two equiprobable states. However, the WT-MCoV complex showed the most conformational variability. Furthermore, to explore where the conformational diversity comes in the complex, the individual RMSD profile of PLpro and hISG15 from each complex was monitored and is shown in Figure B,C. It reveals that the PLpro of each complex remains relatively stable while most of the deviation comes from hISG15.

Figure 2

Distribution of RMSD of (A) complexes; (B) PLpro; and (C) hISG15. Distribution of radius of gyration (RoG) of (D) complexes; (E) PLpro; and (F) hISG15.

The conformational dynamics of hISG15 were explored by comparing Apo-hISG15 with hISG15 from each complex, and the distribution of RMSD is shown in Figure C. It depicts that the Apo-hISG15 is quite dynamic, as seen from the large RMSD value peaks at ∼9 and ∼12.5 Å. In contrast, the highest probable peaks obtained for WT complexes were around 2.5 to 3.0 Å. The SCoV2(MT) and SCoV(MT) showed another likely conformational state at ∼5 Å. However, MCoV(MT) showed the most conformational dynamic states of the hISG15 close to Apo-hISG15, suggesting the more dynamic behavior of MCoV(MT).

In addition, the time-series RMSD of Apo-hISG15 for the entire 4 × 4 μs is shown in Figure S3. It depicts that the Apo form of hISG15 is highly dynamic in nature. The individual N-domain of hISG15 exhibits more deviations than the C-domain, as shown in Figure S3B,C, respectively. Overall, the RMSD trajectory analysis reveals that the Apo-ISG15 is quite dynamic, and it gets stabilized upon interactions with the PLpro’s of CoVs. Moreover, the charge-flip mutation in PLpro also has a significant impact that affects these interactions’ stabilization. Overall, different CoV-PLpro behave differently upon interaction with hISG15, and among them, the MCoV-PLpro-hISG15 complex was found to be the most dynamic system.

Next, we estimated the radius of gyration (RoG or Rg) to measure the structural compactness of each protein–protein complex. The distribution of RoG for each complex, PLpro, and hISG15 is displayed in Figure D–F. In agreement with the RMSD plot, the most significant dynamic states were observed for MCoV compared to others as the distribution is characterized by multiple peaks (Figure D). The RoG of PLpro from each complex shows similar conformational states between the WT and their mutants for each complex, except for MCoV, which exhibited two equiprobable states. However, the major difference in the angular motion was observed for the hISG15, as shown in Figure F. Compared to others, the lower RoG of Apo-hISG15 suggests that, alone hISG15 contributed to low conformational shifts concerning N- and C-domains. In contrast, complexes of PLpro with hISG15 exhibited a large value of RoG ∼18–19 Å, which could be attributed to conformational shifts resulting in changes in the N- and C-domain orientation of hISG15 upon interaction with PLpro during the MD simulations.

We also measured the SASA of all six complexes as well as for Apo-hISG15, and shown in Table S1. The average SASA values varied between 21,721.49 Å2 (SCoV2) and 22,733.31 Å2 (MCoV) for all complexes. The most remarkable change in the SASA value after the mutation was observed in the case of SCoV2 (489.57 Å2). However, in cases of SCoV and MCoV, it decreased by almost 350 Å2. An increase in the SASA value for SCoV2 after mutation indicates the lesser compactness of the structural interface between PLpro and hISG15, which further influences the protein–protein interaction. On the other hand, Apo-hISG15 had a SASA value of 8520.82 Å2, suggesting ∼500–700 Å2 surface area is occupied by PLpro in the process of complexation.

To estimate the dynamical behavior of hISG15 in the free and bound forms, the RMSFs of Cα atom for PLpro and hISG15 were calculated over combined trajectories and shown in Figure BA. In the case of PLpro, a similar pattern of fluctuation was observed for all cases except MCoV and SCoV2(MT). A relatively higher RMSF value was observed in the neighboring region of the mutation site (R166E), especially in the case of SCoV2(MT) (Figure A). The palm’s flexibility is heavily increased in the MCoV’s PLpro, which was subsequently reduced after the point mutation. Compared to all other domains, the thumb domain stays relatively less flexible, which is also observed in other studies.[60,61] Both N- and C-domains of hISG15 in the unbound state showed similar fluctuation, which reduced after the complex formation in both SCoVs and SCoV2, showing a compact interaction compared to other cases.

Figure 3

(A) Root-mean-square-fluctuations (RMSF) of PLpro of each complex system; (B) RMSF of hISG15 of each complex system along with the Apo hISG15.

Interestingly, both domains of hISG15 in the SCoV2 complex show higher rigidity, which is entirely disturbed by the mutation. Comparing all the complexes, hISG15 with SCoV complexes stays in the most stable form. hISG15 in both MCoV and MCoV(MT) shows a higher degree of fluctuation, and also differential flexibility was observed in terms of both N- and C-domains.

Dynamics of hISG15 in the Free and Complex State

The N- and C-domains of hISG15 induce intrinsic flexibility, and hence a critical inspection of its domain motion with respect to each other is needed. The structural stability of both domains was calculated for Apo and complexes by estimating the RMSD distribution. Both the domain’s RMSD was stable throughout our simulation length (Figures S3B,C, S4, and S5). So the dynamics of hISG15 are mostly contributed by the hinge region and interdomain distance. To estimate the two domains’ motion, we calculated the hinge angle for all complexes and in Apo hISG15 and shown in Figure . For calculating the angle, the center of mass of the N-domain, the center of mass of the median amino acids, and the center of mass of the C-domain were used as shown in the red ball and stick model image mentioned above. The free simulation of Apo-hISG15 yields a large conformational basin, including a closed to opened (high energy) structure. However, the closed structure is favorable in the free state (∼70°) (Figure A). Upon binding with PLpro, the global free energy minimum conformation of hISG15 shifted from closed to semiopened or opened state in all wild type and mutant complexes, and the conformational basin shrunk due to the protein–protein interaction. In SCoV2, a mutation in PLpro leads to broadening of the hinge angle, whereas in SCoV(MT), the opposite effect was observed. The global minimum of hISG15 for the complex of SCoV2 was observed around 97° (Figure C) surrounded by a high energy barrier, which got flattened (Global minimum: 111°) in the case of a mutant variant of PLpro (Figure D). On the contrary, SCoV(MT) moves to a semiopened structure (82°) with a higher deviation from its wild type complex (Figure B). In MCoV, both its WT and mutant complexes showed similar trends as SCoV. Thus, the point mutation in the PLpro structure in all three viruses influences the complex formation with hISG15, which changes the two domains’ dynamics.

Figure 4

Potential of mean force of the hinge angle calculated at 300 K. Angle is constructed between the center of mass of N-domain, the center of the hinge region, and the center of mass of C-domain as shown in red dots in the 3D conformation of hISG15. (A) Closed structure of hISG15 in free simulation. (B) and (C) Semiopened structure of hISG15 in complex simulation with SCoV(MT) and SCoV2 PLpro, respectively. (D) Open structure of hISG15 in SCoV2(MT) simulation.

Along with the interdomain movement, we also estimate the two domains’ relative motion with respect to PLpro. For both domains, the center of mass distance between PLpro and the respective domain was calculated. The combined probability distributions are shown in Figure , along with several representative structures of the complexes. The C-domain mostly contributes to the binding of hISG15 with the finger and palm region of PLpro of CoV. N-domain, which is hanging through the hinge region, also makes loose contact with the thumb region. By analyzing the respective distance of each domain of hISG15 with PLpro, one may explore the behavior of hISG15 interaction with PLpro. As depicted by the distribution plot (Figure B), for all cases, the hISG15-domain to PLpro distance was decreased by 5 Å in the C-domain compared to N-domain (Figure A), suggesting strong interactions of PLpro via the C-domain of hISG15. Depending upon the types of PLpro and its mutation, the dynamicity of hISG15 was changed. In SCoV2, two similar populations were observed where both configurations were separated by 5 Å (Figure A). However, in the mutation (SCoV2(MT)), two conformations were observed, but both were shifted closer to PLpro. A similar observation was found in cases of SCoV and SCoV(MT). Contrary to the SCoV2’s bimodal distribution, SCoV and MCoV prefer a unimodal distribution with a higher occupancy in the case of mutation (see Figure C,E). In WTs, the C-domain of hISG15 was comparatively stable in most cases except for SCoV2, where three distinct peaks were observed at 22, 24.5, and 28.3 Å. It indicated a high flexibility level, suggesting that the sampled conformations occupied both close and far states from the PLpro (Figure H). However, the mutant SCoV2(MT) showed that the C-domain moved far from the PLpro with a broader distribution minimum of around 27.5 Å. In the case of SCoV, steady-state interactions were observed between both domains and PLpro, which lie parallel with the PLpro (Figure G). Its mutant variant (SCoV(MT)) exhibited broader distribution, suggesting that it moved closer to PLpro. Almost close contact was also observed with the WT-MCoV complex (Figure F) and even in its mutant MCoV complex. It suggests that the dynamic rearrangement of the respective domain of hISG15 occurs upon interaction with the PLpro of CoVs.

Figure 5

(A) Distribution of the center of the mass distance between PLpro and (A) N-domain of hISG15 and (B) C-domain of hISG15. Structures of several distinct peaks are shown in ribbon representation. (C) MCoV(MT); (D) SCoV2; (E) SCoV(MT); (F) MCoV; (G) SCoV; (H) SCoV2.

PCA of hISG15 in Apo and Complexes

The trajectories of hISG15 from its free simulation and complexes were used in the PCA by taking the Cα atom into account. The eigenvalues corresponding to eigenvectors were calculated and displayed in decreasing order and shown in the Supporting Information (Figure S6). The first few eigenvectors described the collective motion of significant fluctuations, and movements were not the same for all six complexes and Apo-hISG15. The first three PCs contained 84.30, 90.07, and 76.09% of the overall motion for the complex of MCoV, SCoV, and SCoV2. In contrast, the motions were increased in mutant complexes as calculated to be 92.45, 92.95, and 88.04%, respectively. In the Apo simulation, 86.62% of the fluctuation was estimated by the first three PCs. Subsequently, we have constructed a two-dimensional (2D) free energy landscape (FEL) for each system, taking the first two PCs (PC1, PC2) as collective variables, and shown in Supporting Information Figure S7. In all cases, several minima were observed, which indicates the intrinsic flexibility of hISG15. However, the sampling space was reduced in the complex compared to Apo. Still, several low energetic structures were observed for these complexes. In Apo hISG15, three distinct minima separated by almost 3 kcal/mol were observed (Figure S7A). In SCoV2, upon mutation in PLpro, the sampling space increased, and the global conformation position also shifted (Figure D & G). Only in the case of MCoV(MT), we obtained several minima sets in two groups separated by a high energy barrier (more than 4 kcal/mol). It leads to a highly variable structure comparable to the Apo state.

The N-domain and C-domain’s internal structural behavior was almost stable as we already discussed (Figures S3B,C, S4, and S5), so the Cartesian coordinate-based PCA analysis will not be able to capture the entire dynamics. The linker loop between both domains is mostly responsible for the dynamics of hISG15. Hence, we performed dPCA of the hinge region of hISG15 as shown in Figure . This enables us to investigate the dynamics behind the hinge region. The sampled conformational space depicted for Apo hISG15 differs from those of complexes. Also, the difference can be noticed between WT and the corresponding mutant complex. To visualize the structural changes, the k-mean clustering was done corresponding to the hinge region and shown along with the dPCA plot (Figure ). Ramachandran plots were also generated for all significant structures obtained from the dPCA analysis (see Supporting Information Figure S8) to support the structural validity of the same. In the case of hISG15, Figure shows four prominent structural changes in the hinge loop from a bent loop (47.5%) to a linear structure with a lesser probability (11.2%). Compared with the SCoV2(MT) complex dPCA plot, SCoV2 (Figure D) showed higher conformational space variations due to small bend or twist in the loop. In the mutant complex (SCoV2(MT)), hISG15 exhibited two distinct states, mostly by an angular bend near the domain region. This loop structure behavior agrees with the angular distribution of hinge angle (Figure ), where SCoV2 and SCoV2(MT) complexes showed steep and flattened free energy well, respectively. The SCoV complex showed that the two minima conformations were separated by a high energy barrier. At the same time, its mutant form (SCoV(MT)) exhibited that their two conformational states were approachable. A higher degree of variation was observed in the MCoV complex, as agreed with the hinge angle distribution exhibiting broader free energy well. Comparably, MCoV(MT) showed distinct conformations among all other complexes, such as the secondary structure formation in the second-highest conformational cluster (26.3%). Overall, as observed for all complexes, hISG15 exhibits multiple minima, separated by a high energy barrier implying its dynamicity after the complex formation. So, the differential flexibility of hISG15 can be a great measure to estimate the interaction profile and the mechanism of the immune response against a different strain of CoV pathogenesis.

Figure 6

2D FEL generating using the dPCA of the hinge region of hISG15 of each complex obtained from the MD simulations. (A) Apo; (B) MCoV;(C) SCoV; (D) SCoV2; (E) MCoV(MT); (F) SCoV(MT); (G) SCoV2(MT).

Binding Free Energy between PLpro and hISG15

To further evaluate the relative binding free energy between PLpro and hISG15 for MCoV, SCoV, and SCoV2, ΔGbind was calculated based on the MM/PBSA scheme. Moreover, the value of their binding free energy of WT was compared with their respective mutants (R167E). As shown in Table , the various components of the total binding free energy indicate that the intermolecular electrostatic interactions (ΔEele), van der Waals interactions (ΔEvdW), and nonpolar solvation free energy (ΔGnp) favor the binding of PLpro with hISG15. In contrast, the polar solvation free energy (ΔGpol) disfavors the complexation.

Table 1

Binding Energy and Its Various Components for the PLpro of MCoV, SCoV, and SCoV2 with the hISG15 Calculated Using the MM/PBSA Schemea

complex	ΔEvdW	ΔEelec	ΔGpol	ΔGnp	ΔEMMb	ΔGsolvc	ΔGTotald
MCoV	–98.19 (0.11)	–206.62 (0.51)	247.62 (0.46)	–12.81 (0.01)	–304.81 (0.55)	234.81 (0.45)	–70.01 (0.15)
MCoV(MT)	–113.05 (0.11)	–417.66 (0.83)	454.26 (0.72)	–14.44 (0.01)	–530.71 (0.90)	501.98 (0.80)	–90.97 (0.15)
SCoV	–90.56 (0.14)	–229.38 (0.61)	256.44 (0.59)	–12.43 (0.01)	–319.86 (0.65)	244.01 (0.58)	–75.85 (0.17)
SCoV(MT)	–130.45 (0.10)	–142.29 (0.62)	221.68 (0.60)	–14.17 (0.01)	–272.74 (0.64)	207.51 (0.59)	–65.22 (0.19)
SCoV2	–101.35 (0.16)	–231.60 (0.78)	263.92 (0.68)	–12.09 (0.02)	–332.95 (0.82)	251.83 (0.67)	–81.11 (0.19)
SCoV2(MT)	–78.44 (0.18)	–415.64 (0.65)	432.75 (0.66)	–10.04 (0.02)	–494.08 (0.72)	422.70 (0.65)	–71.38 (0.13)

The standard error of means is given in parentheses.

ΔEvdW + ΔEelec.

ΔGnp + ΔGpol.

ΔEvdW + ΔEelec + ΔGnp + ΔGpol.

The SCoV2-PLpro binds more strongly to hISG15 as it exhibited relatively stronger affinity (ΔGbind = −81.1 kcal/mol) than other WT complexes. The predicted ΔGbind for SCoV and MCoV was −75.85 and −70.01 kcal/mol, respectively. It is evident from Table that both ΔEvdW and ΔEelec are more favorable in the case of SCoV2 compared to MCoV or SCoV. This leads to stronger binding between SCoV2-PLpro and hISG15 than the other two variants.

Further, the obtained ΔEMM revealed that for each complex, the electrostatic attraction is the main force inducing the molecular complexation between PLpro and hISG15. ΔEelec ranges from −206.62 to −231.60 kcal/mol while ΔEvdW varies between −90.56 and −101.35 kcal/mol for all WT complexes. It suggests that ΔEelec is ∼2–3 times stronger than ΔEvdW in cases of WT complexes. For each mutant complex, ΔEelec is ∼1–5 times more favorable than ΔEvdW. Notably, the high electrostatic contribution toward the complexation arises from the positively charged arginine (R167) of hISG15.

Next, we investigated the effect of mutation on the PLpro-hISG15 binding. In cases of SCoV(MT) (−65.22 kcal/mol) and SCoV2(MT) (−71.38 kcal/mol), the binding free energy ΔGbind decreases compared to the respective WT. In contrast, the mutation enhances the complexation between PLpro and hISG15 for MCoV(MT) (ΔGbind = −90.97 kcal/mol). It is to be noted that in the case of SCoV2(MT), ΔEelec (−415.64 kcal/mol) increases significantly, while ΔEvdW (−78.44 kcal/mol) decreases compared to WT (ΔEelec = −231.60 kcal/mol, ΔEvdW = −101.35 kcal/mol). However, the net electrostatic energy (ΔEelec + ΔGpol) is less unfavorable to the complexation in SCoV2(MT) than SCoV2. Consequently, the mutation-induced reduced binding affinity in SCoV2(MT) arises primarily due to weaker van der Waals interaction.

On the other hand, in the case of SCoV(MT), the opposite trend is observed where ΔEvdW (−130.45 kcal/mol) increases while ΔEelec (−142.29 kcal/mol) decreases compared to WT (ΔEelec = −229.38 kcal/mol, ΔEvdW = −90.56 kcal/mol). Further, the net electrostatic energy (ΔEelec + ΔGpol) is more unfavorable to SCoV(MT) complexation than SCoV2. This again leads to an overall reduced affinity for SCoV(MT) compared to SCoV, and the weaker electrostatic interaction is mainly responsible for the reduced binding free energy.

Finally, in the case of MCoV(MT), both ΔEelec and ΔEvdW contribute more favorably than WT, leading to an increased binding affinity between PLpro and hISG15. Thus, it justifies that mutation significantly impaired some interactions between PLpro and hISG15, affecting the overall binding affinity. This finding is in good agreement with the reported binding of SCoV-PLpro and hISG15.[32] In addition, our predicted results also agree with different host substrate preferences where it was reported that SCoV preferentially cleaves the ubiquitin chains, whereas SCoV2 predominantly targets ISG15, as seen with mISG15.[62]

Hot-Spot Residues of the PLpro and hISG15

To gain insight into the molecular basis of the specificity of SCoV2 with hISG15, we investigated the crucial residues involved in the binding (ΔGbindresidue) using the MM/GBSA method. The per-residue decomposition of binding energy was estimated and is shown in Figure . Here, we highlighted those residues having free energy values of ≥1.5 kcal/mol and listed them in Table S2. The positive and negative values of ΔGbindresidue represent the stabilization and destabilization, respectively. The more the negative value, the more crucial the residue is for the complexation.

Figure 7

Per-residual contribution to free energy of PLpro and hISG15 complexes of the WT and mutants of MCoV, SCoV, and SCoV2, respectively, obtained via MM/GBSA analysis.

Notably, the number of critical residues from PLpro of SCoV2 was higher compared to other CoVs. It includes Tyr268, Phe69, Pro223, Thr225, Met208, Asn267, Thr74, Thr75, and Asn128, which showed primary interaction with the N-domain residues of hISG15 (mainly with Arg57 of hISG15). In the case of SCoV, the key residues include Val226, Leu76, Asp230, Arg167, Pro224, and Tyr269. This result agrees with the experimental findings of SCoV interactions with the C-domain hISG15 (chISG15).[1] Moreover, it was reported that SCoV showed preferential ubiquitin activity from its Leu76 mediated hydrophobic interaction with Ile44 on ubiquitin. Similarly, we observed from the MD simulations that Leu76 of SCoV favorably interacts with hISG15. This indicates that the hydrophobic residue at this site is crucial in SCoV for its ubiquitin and hISG15 binding.[62] In contrast, the corresponding Thr75 of SCoV2 exhibited comparatively less contribution toward binding with hISG15. In addition, it was earlier shown that the PLpro of SCoV and SCoV2 displayed a similar binding mode with the C-domain mISG15.[32] However, our simulations result predicted the difference in SCoV and SCoV2-PLpro binding behavior with the full-length hISG15. It was observed that the SCoV-PLpro binds more strongly with the C-domain of hISG15 (Arg153 and Arg87), while SCoV2-PLpro binds more favorably with the N-domain (Arg57).

Similarly, the major contributing residues for MCoV are also listed in Table S2. However, interesting results were observed for the mutant SCoV2, where the significantly contributing residues of PLpro lose their interactions with hISG15, leading to an overall decrease in the number of contributing residues compared to SCoV and MCoV. Instead, it is noted that in SCoV2(MT), Arg57 from hISG15 contributes (−11.48 kcal/mol) more than its WT (−6.90 kcal/mol). This illustrates the difference in the interaction profile between PLpro and hISG15 upon mutation. Also, as reported by Daczkowski et al.,[32] the charge-flip mutation of R167E in SCoV-PLpro showed 20 times less efficiency than WT in hydrolyzing ISG15. Similarly, we found that Arg167 contributes to the WT PLpro-hISG15 complexation while replacing the positively charged arginine with negatively charged glutamate abolishes the contribution. It also supported our findings of lower binding affinity observed in SCoV(MT) than SCoV.

Similarly, apart from the closely related family of SCoV and SCoV2, the diverse family of MCoV showed that the charge-flip mutation increases the number of different contributing residues, such as Arg87, Arg92 from hISG15 and Arg82 from the PLpro. It may enhance some more interactions leading to its high affinity than its WT.

PLpro-hISG15 Interactions

Since the electrostatic interaction was the main force behind the PLpro-hISG15 complexation, we investigated the intermolecular hydrogen bond (H-bond) formation using the MD trajectories. The H-bond occupancy (%) is recorded in Table S3 and determines the strength and stability of the interaction. Also, the H-bond interaction between the WT or mutant SCoV2 PLpro and hISG15 is shown in the ball and stick representation to understand the key changes due to the point mutation (Figure S8A,B).

It is found that hISG15 formed H-bonds with several polar and charged residues of CoV-PLpro. Notably, the number of residues forming a strong H-bond between PLpro and hISG15 (>50% occupancy) was higher in SCoV2 than other CoVs. As seen from Table S3, Arg57 of hISG15 forms strong H-bonds with PLpro residues, including Gln174 (∼51%), Asn128 (∼50%), and Asp179 (∼40%). Similarly, other H-bond interactions include Thr74-Ser21 (∼48%) and Thr225-Pro144 (∼35%). In contrast, in the case of SCoV2(MT), the H-bond interaction decreases as seen with Gln174-Arg57 (∼13%), while Arg57 from hISG15 forms another interaction with Glu166 (two H-bonds, ∼52%) and Glu167 (one H-bond, ∼21%) (see Figure ). Similarly, the main H-bond formed in SCoV was between Arg87 of hISG15 and Asp230 (two: H-bond, ∼41%). The other important H-bond was formed between Thr147@hISG15 and Val226@PLpro (∼34%). However, in the case of SCoV(MT), the primary H-bond between Asp230-Arg87 was not persistent (∼16%). Instead, other H-bonds with significant occupancy were formed, such as Met209-Gly128 (∼54%), Gly228-Thr147 (∼27%). In the case of MCoV, the significant H-bond interactions were formed between Gln227-Leu145 (∼53%), Thr147-Val225 (two: ∼50%, 35%), and Thr147-Tyr224 (∼32%). Upon mutation (MCoV(MT)), the Gln227 residue forms the main H-bond interaction with the Ser88 residue of hISG15 (∼42%). Other interactions include Phe128-Ser26 (38%) and Glu178-Ser31 (∼24%). Together, it reveals that the nitrogen atom of Arg57 from hISG15 forms hydrogen bonding with the oxygen of glutamine (Q174) of PLpro in SCoV2, which upon the charge-flip mutation (R166E), switch to form with the negatively charged glutamate (E166).

Figure 8

H-bonds and hydrophobic interactions between the hISG15 and PLpro using the Ligplot (dimplot) for the WT complexes of (A) MCoV; (B) SCoV and (C) SCoV2; respectively. H-bonds are displayed in green dotted lines, and red spikes represent the hydrophobic interactions. Significant H-bonding residues of hISG15 and PLpro are highlighted in green and blue, respectively.

In addition to the hydrogen-bond interaction, the hydrophobic interaction pattern governing the complex formation between hISG15 and PLpro is investigated and drawn using the Ligplot, as shown in Figure and Figure S9. The WT SCoV2-PLpro showed many H-bonds and hydrophobic interactions compared to SCoV and MCoV. Hydrophobic interactions are observed, including residues like Arg153, Phe149, Thr147, Ser146, Gly128, Val59, Ser22, Pro130, Pro59, Met23, and Trp3 from hISG15 and residues Gly209, Pro233, Ala176, Cys224, His175, Ile222, Met208, His73, and Phe69 from PLpro. This indicates a strong tendency of binding of SCoV2-PLpro with hISG15. In contrast, the mutant MCoV showed a large number of H-bonds and hydrophobic interactions compared to others. These results also supported the high binding energy for MCoV(MT) toward hISG15 than WT. Overall, our study provides evidence that SCoV2-PLpro shows preferable binding toward hSIG15 than other CoVs.

Residual Network Analysis

The NAPS server provided a visual inspection of the subnetwork based on various physicochemical properties, such as hydrophilic, hydrophobic, and charged interactions. The different complex nodes between PLpro and hISG15 of all systems were explored and are listed in Table S4. The result depicts the significance of the hydrophilic interaction network in WT-SCoV2 compared to other complexes, such as SCoV and MCoV (Figure ). No hydrophobic or charged nodes contact was detected for WT-SCoV2. Instead, the important hydrophilic interaction networks were found between PLpro and hISG15, which includes (Thr225-Thr147), (Gln174-Gln63), (Asn128-Ser62), (Thr74-Ser24, Ser22), and (Thr75-Ser21, Ser22). However, for SCoV2(MT), a lesser number of hydrophilic network nodes was observed. Additionally, some hydrophobic nodes were also observed, mainly (Leu199-Val52, Ala53) (Figure S9). In contrast, WT-SCoV exhibited mainly hydrophobic contact nodes, while upon R167E mutation in SCoV-PLpro, it showed more hydrophilic and charged interaction networks (see Table S4 and Figures S10 and S11). Similarly, for the WT-MCoV complex, mainly hydrophilic network nodes were prominent, while its mutant form displayed an increase in both hydrophobic and hydrophilic contact nodes compared to charge nodes. In general, it can be concluded from residual network analyses that SCoV2-PLpro interacted heavily with the hydrophilic residues of hISG15 resulting in high binding affinity.

Figure 9

Significant hydrophillic nodes (red) obtained from the network analysis between PLpro (blue) and ISG15 (green) of (A) MCoV; (B) MCoV(MT); (C) SCoV; (D) SCoV(MT); (E) SCoV2 and (F) SCoV2 (MT); respectively. Important connections are shown in yellow connective lines.

Conclusions

This comprehensive study offers an integrated computational approach to explore the detailed mechanism of all three major CoVs-PLpro’s interaction with the complete human ISG15 for the first time. The PLpro of CoVs forms the replicase complex and antagonizes the innate immune responses via deubiquitinating and deISGylating activities, making it a potential target in treating the COVID-19. In the absence of the complete co-crystallized structure of hISG15 and PLpro, we constructed the same, and MD simulations were performed. Because of the bi-lobular construct of hISG15, each domain’s respective orientation and their differential interaction pattern with the conserved CoVs were investigated. By analyzing the dynamics of the free hISG15 protein and bound with PLpro of CoVs, we have provided insights into the differences in substrate preferences, binding mode, and the intermolecular interactions relevant to the function of SCoV2, SCoV, and MCoV. Despite the same sequence, hISG15’s dynamics of both domains were differently influenced upon binding with PLpro of CoVs. Our binding free energy analysis suggests that hISG15 binds more strongly with SCoV2-PLpro compared to other CoVs, and mainly the electrostatic interaction drives the complexation between hISG15-PLpro. This may explain why SCoV2-PLpro preferentially cleaves hISG15, whereas SCoV-PLpro predominantly targets ubiquitin chains. The C-terminal domain of hISG15 mainly governs all the interactions between PLpro and hISG15, while the N-terminal domain interacts weakly with PLpro. It agrees well with the recent experimental studies.[33] These in-depth analyses would fill the critical gap of understanding how PLpro of the CoV family interacts with hISG15, leading to specificity in SCoV2. This study may help in designing drugs targeting PLpro of SARS-CoV-2 that can suppress the COVID-19 infection and promote antiviral immunity.