GRP78: A possible relationship of COVID-19 and the mucormycosis; in silico perspective.

Mucormycosis is a severe fungal infection reported in many cancer survivors, diabetic and immune-suppressed patients during organ transplants. A vast spark in the reported COVID-19 cases is noticed in India during the second wave in May 2021, when Mucormycosis is declared an epidemic. Despite being a rare disease, the mortality rate associated with Mucormycosis is more than 40%. Spore coat proteins (CotH) are essential proteins in many pathogenic bacteria and fungi. CotH3 was reported as the vital protein required for fungal virulence in Mucormycosis. We previously reported the involvement of the host cell-surface receptor GRP78 in SARS-CoV-2 spike recognition. Additionally, GRP78 is known to be the virulence factor during Mucormycosis. Using state-of-the-art structural bioinformatics and molecular modeling tools, we predicted the GRP78 binding site to the Rhizopus delemar CotH3 protein. Our findings pave the way toward rationally designing small molecule inhibitors targeting the GRP78 and its counter proteins in both pathogenic viral (SARS-CoV-2 spike) and fungal (R. delemar CotH3) diseases.

Introduction

Mucormycosis, previously termed zygomycosis, is a severe but rare fungal disease caused by Mucorales order [1]. Rhizopus oryzae is the primary representative organism causing the illness, as observed in ∼70% of mucormycosis cases [2,3]. Mucormycosis is frequently reported after massive natural catastrophes such as tsunami and tornados [4]. Therefore, many patients can be at risk of mucormycosis including patients with uncontrolled diabetes mellitus, in ketoacidosis, undergoing organ or bone marrow transplantation, treated with corticosteroids, have trauma and burns or have malignant hematologic disorder [5,6]. Mucormycosis is a fatal infection associated with a high mortality rate, particularly for patients suffering from diffused diseases [4,7]. As mucormycosis usually occurs in immunocompromised hosts, it has been reported recently as a secondary infection for SARS-CoV-2 patients [8], with a massive surge of mucormycosis characterizing the second wave of COVID-19 disease in India [[8], [9], [10]].

The 78-kDa cell-surface glucose-regulated protein CS-GRP78 (also known as BiP and HSPA5) plays a crucial role in the virulence of Rhizopus fungi, which causes mucormycosis infection. GRP78 is recognized by the fugal spore coat proteins (CotH - mainly the CotH3), which act as fungal ligands GRP78 [11]. It is also important to note that we have recently reported the importance of CS-GRP78 in recognizing the SARS-CoV-2 spike protein [[12], [13], [14], [15]], where these findings were validated experimentally [16]. Moreover, GRP78 has been hypothesized to be responsible for the cross-vaccination reported for human coronaviruses [12,15]. The recognition of SARS-CoV-2 spike of the new variants alpha (UK variant VOC-202012/01), beta (South African 501.V2), and Gamma (Brazilian B.1.1.248) by cell-surface GRP78 is enhanced compared to the wildtype SARS-CoV-2 [17, 18].

Growing evidence is emerging to interconnect COVID-19 infection with mucormycosis [19]. ACE2 is the main entry gate for SRARS-CoV-2, but GRP78 recognizes the spike as an auxiliary route of infection, as proved experimentally lately by Carlos et al. [20]. Additionally, the GRP78 is involved in the translocation of ACE2 to the cell membrane. Therefore, mucormycosis susceptibility is increased in stressed cells, like the case in diabetes mellitus, cancer, and viral infection (such as COVID-19).

Here, we are using state-of-the-art computational bioinformatics tools to understand the link between these two deadly illnesses. Protein modeling followed by molecular dynamics simulation and protein-protein docking has been adopted for the fungal spore coat protein (CotH3). At the same time, its binding affinity and mode of binding against CS-GRP78 are predicted.

Materials and method

Sequence retrieval and alignment

Sequence for the fungal spore coat protein CotH3 (RO3G_11882) was retrieved from the National Center for Biotechnology Information (NCBI) sequence database, while the Pep42 sequence was retrieved from the literature [[21], [22], [23]]. RO3G_11882 is the CotH3 protein sequence of Rhizopus delmar, one of the reported potential fungal agents that causes Mucormycosis in humans [24]. Also, the bacterial (Bacillus cereus) CotH and the Eukaryotic Protein Kinase (EPK) sequences were downloaded from the Protein Data Bank (PDB) database [25] (PDB ID: 5JD9 [26] and 1ATP [27], respectively). Sequence alignments were performed using the Clustal Omega multiple sequence alignment web server [28] of the European Molecular Biology Laboratories-European Bioinformatics Institute (EMBL-EBI). ESPript 3.0 web server was used to visualize the alignment [29].

Structure prediction and the docking of AMP

The fungal CotH3 all atoms 3D model was built using the homology modeling SWISS-MODEL web server of The Center for Molecular Life Sciences, University of Basel [30]. We modeled the part of the CotH3 that is covered by the bacterial homolog (PDB ID: 5JDA). This model is missing 171 residues from the N-terminal and 133 residues from the C-terminal. This region of the CotH3 that we modeled has the whole CotH region of interest in our study (CotH kinase protein) as per the NCBI protein database definition of the protein regions of RO3G_11882. The generated model was validated using the SWISS-MODEL tools for validation and by the Structural Analysis and Verification Server (SAVES) of the University of California, Los Angeles (UCLA) [31]. Three software were used from this server in model validation, including PROCHECK [32], ERRAT [33], and Verify-3D [34]. The active site amino acids for the fungal CotH3 are H135, R208, K266, Y280, and Q386 based on sequence and structural alignment against the bacterial CotH (PDB ID: 5JDA) [35].

In eukaryotes, the activity of the kinase is managed by nucleotide-binding such as Adenosine monophosphate (AMP) and the ATP: AMP ratio [36]. Nucleotide-binding is reported to be responsible for the kinase activity of human and bacterial kinases as well [35,37]. In this study, we docked AMP to the active site of the fugal CotH3 to test its binding activity against the host cell receptor GRP78. AutoDock Vina software was used to perform AMP docking against the CotH3 model, where we docked it to the protein active site [38]. A flexible ligand in a flexible active site docking protocol was adopted, where the formed complexes were ranked by their estimated Vina scoring function [38]. The docking search space used a grid box of size 60 Å × 90 Å × 50 Å, centered at (35.4, 60.5, 43.3) to cover the active site of CotH3. After that, two magnesium ions are added to the model using the coordinates of the bacterial CotH structure (PDB ID:5JDA) after superposition with the fungal CotH3 model.

Molecular dynamics simulation (MDS)

The fungal CotH3 with AMP and two Mg+2 system was solvated in a TIP3P water box at pH 7 [39]. The total charge of the system was maintained at zero by adding 41 sodium ions and 30 chlorine ions to mimic the physiological NaCl concentration of 0.154 mol/L. A total number of 10,901 water molecules in a box of size 70.47 Å × 78.27 Å × 70.92 Å were added, making a whole simulated system of 35,640 atoms. Water molecules in the system were then minimized for 10000 steps using a conjugate gradient algorithm [40]. The entire system (CotH3, AMP, Mg+2, and water) was then minimized using the same algorithm for another 10000 steps. After that, the constrained system (protein atoms fixed) was heated up to 310 K (physiological temperature), followed by an equilibration MDS run for 100 ps at NVT ensemble (constant number of atoms, volume, and temperature). Subsequently, an equilibration run (1 ns) for the whole system, without constraints, was performed. Finally, a production run of 120 ns was performed using the same conditions of the equilibration run for the fungal CotH3 system. Root Mean Square Deviation (RMSD) and Radius of Gyration (RoG) were plotted to show the equilibration of the system versus the simulation time. Hydrogen bonds (H-bond) formed between AMP and the protein were recorded during the 120 ns MDS run.

Protein-protein docking experiment

After the production run, the MD trajectories were clustered into groups of similar conformations using the Chimera software [41]. The clustering was performed for the trajectories after 60 ns (50% of the trajectories) of the MDS to ensure system equilibration. Ten different conformations representing ten clusters were used to test the binding affinity of the GRP78 (PDB ID: 5E84, Chain A) against the CotH3 model using HADDOCK software 2.4 [42]. Complexes of GRP78-CotH3 were then analyzed using the Protein-Ligand Interaction Profiler (PLIP) webserver to check the formed interactions then tabulated [43]. The PyMOL software was utilized to generate the 3D figures presented in this manuscript [44,45]. After that, a 50 ns MDS run was performed for the GRP78-CotH3 complex using the same protocol to check the dynamics of the formed complex.

Results and discussion

We are dealing with a life-and-death situation by combining a fatal fungus with a widely spread viral infection. This raises an important question as what could be the interconnection between the two diseases? We hypothesize that the stress chaperon protein (GRP78) and spore coating protein H3 (CotH3) provide this link for more than one piece of evidence, as we will see in the following sections. In the current study, we don't say that direct host-virus-fungus links persist; instead, the host cell protein GRP78 has the dual capability to work as an internalization gate for fungal and viral infection elements (CotH3 and Spike, respectively).

Multiple sequence and structural alignment

No structures are found in the protein data bank database for the fungal spore coating protein H; hence we have to predict the 3D structure computationally. But first, Do the available bacterial spore coat protein structures can help!. The structures for both eukaryotic protein kinase and bacterial CotH are available in the protein data bank. Sequence alignment of the eukaryotic protein kinase (PDB ID: 1ATP), the fungal CotH3 (RO3G_11882), and the bacterial CotH (PDB ID: 5JD9) reveals interesting results (see Fig. 1 ).

Fig. 1

Multiple sequence alignment of eukaryotic protein kinase (EPK) (PDB: 1ATP), Fungal CotH3, and Bacterial CotH (PDB: 5JD9).

A number of features in fungal CotH3 resemble those in both the eukaryotic protein kinase (EPK) and prokaryotic kinases. For example, the Glycine-rich region and the APE motif characterize EPK in the fungal CotH3 but are absent in the bacterial CotH [46]. Furthermore, sequence identities between fungal CotH3 and EPK, & fungal CotH3 and bacterial CotH are 17.71% and 19.78%, respectively, while it is only 10.0% between the EPK and bacterial CotH. Based on this multiple sequence alignment, we propose that the fungal CotH3 can function as a kinase-like EPK and the bacterial CotH.

CotH3 model construction

Fungal CotH3 all atoms 3D model (296 amino acids) was constructed using the homology modeling web server SWISS-MODEL (7). The only suitable template with good coverage for the fungal spore coat protein was the Bacillus cereus CotH, despite its low sequence identity (19.93%). The best model constructed for fungal CotH3 was valid based on the results of the Ramachandran plot (98.5% in the preferred region, 1.5% in the allowed region with no outliers) (see Fig. 2 A), Verify-3D (82.15% of the residues had an averaged 3D-1D score greater than 0.2), and ERRAT (overall quality factor is 80.9%).

Fig. 2

A) The Ramachandran plot of the predicted model of CotH3 of Rhizopus delmar. B) The predicted average binding affinity of AMP to the fungal CotH3 (green) and the bacterial CotH (red) apo structure (PDB ID: 5JDA).

For CotH3 to be active as kinase AMP should be present in its position, we docked it using AutoDock Vina software. Fig. 2B shows the average binding energy of AMP to the fungal spore coat protein model (green column). AMP can bind to CotH3 with a binding affinity of (−7.6 to −8.6 kcal/mol). Besides, AMP is docked into the solved structure for the apo form of the bacterial CotH (PDB ID: 5JDA) with a predicted binding affinity of −7.9 kcal/mol (red column) the same range as that for the fungal spore coat protein.

Based on these results, we suggest that the fungal spore coat protein can tightly bind to AMP and may function like the bacterial CotH protein as a kinase, but yet to be validated experimentally.

Molecular dynamics simulation

The fungal CotH3 model was minimized and equilibrated with molecular dynamics simulation for a period of 120 ns. This experiment was conducted to be sure of the different conformations CotH3 will take during this time interval. As shown in Fig. 3 , the system was equilibrated after about 60 ns.as reflected from the Root Mean Square Deviation (RMSD) values (Fig. 2A), where the RMSD (blue line) is stabilized at about 6.5 Å. Also, based on the radius of gyration (RoG) values (red line), the system is stable during the entire period of the MDS run, with the radius of gyration fluctuates between 20 and 22 Å.

Fig. 3

The MDS analysis of the CotH3-AMP system during 120 ns. (A) The Root Mean Square Deviation (RMSD) in blue and the Radius of Gyration (RoG) in red, versus the simulation time (in nanoseconds). (B) H-bond occupancy versus time in ns. (C) The per-residue Root Mean Square Fluctuations (RMSF). Protein conformation is represented in colored cartoon representation according to the coloring scheme on the up-right corner.

Hydrogen bonds formed between the AMP and CotH3 active site pocket amino acids were tracked during the entire period of the MDS run, as shown in Fig. 3B. The percent occupancy of the hydrogen bonds was calculated for every 5 ns. We focused on the H-bonds formed between AMP and G179, K181, N184, N368 & D387. As shown in Fig. 3B, the percent occupancy of the AMP-K181 H-bond was maintained during the entire MDS)ranged from 100 to 200% (. Both AMP-N368 and AMP-D387 H-bonds show 50–100% occupancies during the first 100 ns of the MDS. After that, it was dropped and was compensated by AMP-G179 and AMP-N184 H-bonds. So, the AMP was tightly interacting with the CotH3 by at least 3 H-bonds during the simulation time. This supported the docking study and proved the ability of AMP to be settled in the active site pocket of CotH3.

Fig. 3C shows the per-residue Root Mean Square Fluctuations (RMSF) for the CotH3 system during the simulation time. After the simulation, the 3D structure of the CotH3 model is represented at the top of the figure with colored cartoons. In addition to the protein termini, six regions show high flexibility characterized by RMSF values greater than 3 Å, including the S194-G200 (blue), I243–F248 (cyan), G282-S284 (magenta), G293-A312 (red), L396-N402 (orange), and R415-G418 (gray) regions. The defined kinase activity residues R208, E219, Y281, D360, T367, and D387 for CotH3 (shown on the RMSF curve by location marks) are found to be rigid during the simulation period. The orange region (L396-N402) is suggested to bear the GRP78 binding motif based on the sequence similarity with the peptide Pep42 that was previously reported to selectively bind to GRP78 over cancer cells [22,23,47,48]. Noticeably, the GRP78 recognition site on the spike protein of SARS-CoV-2 (C480–C488) exhibits a similar pattern of elevated RMSF. Additionally, the new strain spikes of SARS-CoV-2 (alpha, beta, and gamma strains) show higher RMSF values for the C480–C488 region than the WT spike. A correlation between the elevated RMSF of this region and the predicted spike recognition by GRP78 is reported as well [17,18]. So, we suggest the CotH3 L396-N402 region to be the recognition element for CS-GRP78 on epithelial cells. Now we will check for the GRP78-CotH3 binding at the predicted sites using the protein-protein docking protocol.

The predicted CotH3 binding site to human GRP78

As reported earlier, the fungal CotH1, CotH2, and CotH3 form direct contact with the human receptor GRP78 found on the cell membrane of the endothelial cells with priority for the binding were for CotH3 [11]. Therefore, we utilized HADDOCK 2.4 webserver to simulate the binding of the major virulence factor for Mucormycosis (CotH3) to the host cell-surface GRP78. The ten representative conformations of the CotH3 after the MDS run were docked to the GRP78 structure. HADDOCK 2.4 utilizes solvated docking to simulate protein-protein interaction using combined information sources such as bioinformatics, NMR, and mass spectrometry to drive the docking. In addition, it includes the depreciation of the protein-protein interface using molecular dynamics [42]. Fig. 4 A shows the average docking scores (HADDOCK scores) for the docking of GRP78 substrate-binding domain β (SBD) to the CotH3 (L396-N402 region) (red column). Additionally, the average docking scores for the wildtype (WT), alpha, and beta & gamma variants of the SARS-CoV-2 spike against GRP78 SBD are shown for comparison (blue, cyan, and purple columns, respectively).

Fig. 4

(A) The binding affinity of the CotH3 (blue column) and Spikes of SARS-CoV-2 (blue, cyan, and purple columns) against GRP78 calculated using HADDOCK software. (B) The docked complexes are superimposed on each other (left-hand side), and one of the formed complexes (right-hand side) shown in cartoon representation (top) and surface presentation (below) with 90° rotation on the x-axis. C) The formed interactions in (B) represented by PyMOL after the PLIP webserver run. The different interactions are depicted as per legend at the bottom of the figure.

The average HADDOCK score for the CotH3 against GRP78 (−76.48 ± 8.0) is in the same range as the different variants of the spike against GRP78 (−85.8 ± 9.8). The formed interactions upon docking the CotH3 against GRP78 are listed in Table 1 . The main types of interactions that stabilize the complexes are the formation of H-bonds (8 ± 2.5) and the hydrophobic interactions (5.7 ± 1.8). Additionally, salt bridges are formed in some complexes (1.1 ± 0.9). The residues from the GRP78 that engaged in the interaction with CotH3, ranked by the number of formed interactions, are V429 (15), S452 (15 H-bonds), R488 (11), V453 (8), T428 (7), T452 (7 H-bonds), V490 (7), T428 (6), T434 (6 H-bonds), K460 (6), Q449 (5 H-bonds), and Q492 (5).

Table 1

The interactions formed between the CotH3 of R. delemar and cell-surface GRP78 SBDβ upon docking with HADDOCK.

CotH3 conformation	HADDOCK score	H-bonding			Hydrophobic interaction			Salt bridges
		number	Residues from the CotH3	Residues from GRP78	number	Residues from the CotH3	Residues from GRP78	number	Residues from the CotH3	Residues from GRP78
1	−66.1±4.9	7	E400D406N402P399E400P465N461	E427V429S452T458K460G489Q492	3	D406L457L464	V432I450V490	1	D466	R488
2	−82.5±3.7	5	Q246N245N395E400D411	D350D350S452T458G489	4	A397N402P399E400	T428V429I459K460	1	E400	K460
3	−73.1±10.6	7	N245P399E468A397D406E400K410	D350V429K453S452S452T458R488	5	L464L396E400A398E468	V432I450V453V453P467	2	E468D411	K435R488
4	−73.2±8.7	5	E400H460D406A398N402	Q449Q449R488Q492Q492	8	N461A458E400A398L457T453L405L405	T428V429I450I450F451V490V490Q492	3	E468E400D411	K435K447R488
5	−91.1±5.5	9	D406N402E400E400E400E400Q216A185G247	G431T434T434S448Q449Q449S452R488G513	9	A398A398N395A397L405A398L396F248E211	I426T428V429V429V432F451V453R488V490	1	E400	K447
6	−70.4±8.9	11	N402N402D406N402E400E400E400F242R240R240L464	T428V429G431T434S452S452V453R488V490V490Q492	5	A397D406P465E400F242	V429V432I450F451R488
7	−78.2±9.6	8	N402N402G401E400E400E400P214K186	T428V429T434T434T434Q449S452V490	4	D406T453A398A185	V429V432F451R488	1	E400	K447
8	−67.0±5.8	6	A398N395N395V394N395D411	V429S452S452S452V453G454	6	L396L396L405A397A398Y416	I426T428V429V429V432P487
9	−87.9±11.2	9	H460H460L396V394D406E400E400K410D411	E427G430S452S452G454T458T458A486G489	7	P399A397L396L396A403A398P399	E427V429F451F451V453V457K460	2	E400D411	K460R488
10	−75.3±7.1	13	N244L396N395E400E400E400N461N461H460H460N244Q246N244	K123G430G430S452S452V453G454T456T458T458E532E532R536	6	F248A398E400E400Q246N245	K123V429F451V453T530E532

Fig. 4B shows the docking complexes formed after docking the GRP78 structure (green cartoons) against the ten different cluster representatives of CotH3 (cyan cartoons). The region F392–V407 of the CotH3 (GRP78 binding region) is shown in the red cartoon. This region is surface accessible and has high hydrophobicity index of 0.663 (Kyte & Doolittle) [49]. This agrees with previous reports that GRP78 can catch unfolded hydrophobic patches on misfolded proteins and viral proteins [12,[50], [51], [52]].

The left-hand side of the figure shows the superposition of the ten complexes, while the right-hand side shows the complex formed between GRP78 and the CotH3 conformation at 82 ns. The surface representation at the bottom view of Fig. 4B shows the binding region of the CotH3 docked into the GRP78 substrate-binding domain β. The PLIP webserver was utilized to analyze the docking complexes (Table 1 and Fig. 4C).

H-bonds (blue lines) are the primary formed interaction types, followed by the hydrophobic interactions (dashed-gray lines) and few salt bridges (dashed-yellow lines with two yellow balls).

We perform 50 ns MDS run on one of the formed complexes (GRP78-CotH3) utilizing the same protocol. Fig. 5 A shows the RMSD in Å (blue line), RoG in Å (orange line), and SASA in Å2 (gray line) for the GRP78-CotH3 complex versus the simulation time in ns. As reflected from the RMSD curve, the system was equilibrated during the first ten ns of the simulation with an equilibrium RMSD value of 8 Å. The RoG and SASA of the GRP78-CotH3 complex indicate system stability during the simulation, with average values of 43 Å and 47000 Å2 for RoG and SASA, respectively. In addition, the total number of H-bonds in the system was stable around 1470, which also indicates system stability. The per-residue RMSF of the GRP78-CotH3 system (red line) was plotted in Fig. 5B along with GRP78 alone (green line) and CotH3 alone (blue line) [53]. The interacting regions in both proteins that we predicted are enlarged for clarification. Overall the RMSF of the GRP78-CotH3 complex (red) is slightly higher than the single protein RMSFs (green and blue). Noticeably, the predicted binding site of the CotH3 to GRP78 (L396-N402) has lower RMSF (red line) compared to the CotH3 RMSF (blue line) as clarified in the top-right enlarged panel. This reflects the stabilization exerted on this loop (orange loop in Fig. 3C) upon binding to GRP78. For GRP78, the RMSF of the complex is close to the free protein as shown from the left-enlarged panel in Fig. 5B.

Fig. 5

Molecular Dynamics Simulation (50 ns) of the GRP78-CotH3 complex after docking. (A) The Root Mean Square Deviation (RMSD) in blue, the Radius of Gyration (RoG) in orange, and the Surface Accessible Surface Area (SASA) in gray, versus the simulation time in nanoseconds. (B) The per-residue Root Mean Square Fluctuations (RMSF) of the GRP78-CotH3 complex (red), GRP78 (green) and CotH3 (blue). The interacting regions in both protein are enlarged for clarification.

In summary, our results don't contradict the ACE2 role in viral (SARS-CoV-2 spike) recognition and entry. Instead, GRP78 over the stressed cells has been proven to be an auxiliary entry element for SARS-CoV-2. Additionally, the membrane expression of ACE2 is elevated only in the presence of GRP78, so its role involves the translocation of ACE2 to the plasma membrane [20]. GRP78 is a stress response inside the cell, so the viral infection will elevate the level of GRP78 expression that translocate more ACE2 and CS-GRP78 to the membrane of the host-cell, increasing the susceptibility of SARS-CoV-2 spike and fungal CotH3 recognition and infection propagation.

Conclusion

Fungal CotH3 is an essential factor for Mucormycosis virulence. In this study, we reported for the first time the kinase activity of CotH3. Besides, we predicted the binding site on CotH3 against the human cell-surface GRP78 protein that is overexpressed on the membrane of endothelial cells upon cellular stress. CS-GRP78 is one of the routes for SARS-CoV-2 recognition and entry, in addition to its master role in recognizing CotH3 and the internalization of the Rhizopus fungal species. Therefore, inhibition of the CotH3-GRP78 binding is a key for suppressing the virulence of Mucormycosis. At the same time, anti–CS–GRP78 may be suitable to reduce the virulence of Mucormycosis-COVID-19 coinfection.

Funding

This paper is partly based upon work supported by Science, Technology & Innovation Funding Authority (STDF) under grant number: 44575.

Declaration of competing interest

All the authors declare that there is no competing interest in this work.