Looking into mucormycosis coinfections in COVID-19 patients using computational analysis

1.   Introduction

In early May 2021, the rare disease mucormycosis has been declared an epidemic and a life-threatening infection in India and other South Asian countries. Mucormycosis, also known as zygomycosis and phycomycosis, is a rare and uncommon infection caused by a group of fungi named mucormycetes that are found to be residing in the atmosphere and primarily affect individuals with comorbidities or immunocompromised health problems. The inhalation of fungus spores is the common mode of contamination, however, it can happen as a consequence of a bodily cut or a burn [1],[2]. There are different types of mucormycosis, including 1) Rhinocerebral mucormycosis which affects the sinus and can spread to the brain, commonly occurring in uncontrolled diabetic patients and people who had kidney transplant [3]–[5]; 2) Pulmonary mucormycosis which is commonly observed in cancer patients and people who had organ or stem cell transplant [6]; 3) Gastrointestinal mucormycosis which is commonly observed in young children and premature born babies who might have surgery or received medications that lowered their germ defense ability [7],[8], 4) Cutaneous mucormycosis which occurs when the fungi break through the skin due to cuts, burns, surgery or skin trauma that can be observed in people who do not have compromised immune system [9]; and 5) Disseminated mucormycosis which spreads through the bloodstream and mostly affects the brain and other organs including the heart and spleen [10].

India experienced the second wave of coronavirus 2019 (COVID-19) infections nearly a year after the pandemic was announced. The dominant strain was then named the Delta variant and classified as a strain of concern due to its increased transmissibility and disease severity, according to the Centers for Disease Control and Prevention (CDC) variants classification [11]. Virus strains that scientists consider are more transmissible or capable of producing more severe illnesses are classified as variants of concern (VOC), including B.1.1.7 lineages (Alpha variants), B.1.351 lineages (Beta variants), P.1 lineages (Gamma variants), B.1.427 and B.1.429 lineages (Epsilon variants), and B.1.617.2 and AY lineages described as Delta variants [12].

Mucormycosis symptoms vary depending on where the fungus is growing in the body and some of the mucormycosis symptoms are similar to those symptoms reported after COVID-19 infection. Symptoms of rhinocerebral mucormycosis can include headache, sinus or nasal congestion, and fever; symptoms of pulmonary mucormycosis may include fever, cough, chest pain, and shortness of breath; while symptoms of gastrointestinal mucormycosis may include abdominal pain, nausea, and vomiting. Mucormycosis is more common in people with diabetes and ketoacidosis, as well as those receiving high dosages of corticosteroids after COVID-19 infection. The combination of various clinical data and the isolation of the fungus from clinical samples in culture is required for the probable diagnosis of mucormycosis [1],[13],[14].

Several risk factors have been reported that would increase the mucormycosis epidemic during the COVID-19 pandemic. Diabetes mellitus, chronic hypertension, cardiovascular diseases, and renal diseases are among the most prominent risk factors linked to an elevated mucormycosis infection rate among COVID-19 patients in India, South Asia, the United States of America, Egypt, Iran, Brazil, Chili [15]–[18]. Patients with a past medical history of one or more comorbidities and those who have recovered from COVID-19 infection are more susceptible to mucormycosis infection. In addition, it has been noted that unwarranted medication used for treating severe COVID-19 symptoms might have increased patients' likelihood of mucormycosis infection, where steroids are the most frequently prescribed medication for COVID-19, followed by Remdesivir, antibiotics, and Tocilizumab [16],[19]. The infection rates have been noted to vary among different regions where high infection rates in the United Kingdom, France, Italy, Austria, and Mexico, is reported among COVID-19 patients with organ transplants and immunocompromised people [20]–[23].

Mucormycosis disease incidence and infection rates have been increasing, particularly during and after infection with the virus causing the COVID-19, and have been observed in India during its second wave of infection [24],[25]. Patients with severe COVID-19 symptoms, particularly those admitted to the hospital and into the intensive care unit, are more prone to develop this fungal infection, which has been associated with serious illness and death [22]. COVID-19 specific treatments, such as receiving high-dosage corticosteroids that had been used to treat severe COVID-19 cases, are more likely to increase the risk of mucormycosis in COVID-19 patients [17],[26].

Mucormycosis is not a new form of fungus, it is also known as the black fungus, and it is caused by species called Mucorales with Rhizopus Oryzae being known the most. Up to 70% of all cases of mucormycosis are caused by this fungus [27]–[29]. Mucormycosis infection is more likely in individuals who have a weak immune system. Patients with comorbidities like diabetes are more susceptible to contract mucormycosis, as are those taking steroids to treat severe COVID-19 symptoms. Mucormycosis, on the other hand, is not a contagious disease [30]. Intravenous catheterization and the use of broad-spectrum antibiotics are risk factors for getting mucormycosis. Mucormycosis can also develop as a result of surgical procedures, hyperalimentation, or malnutrition [31].

Global epidemiological studies of mucormycosis are reasonable in assessing the disease pattern and the infection incidence among people who are at high risk in different countries. Mucormycosis has diverse causal agent factors depending on the geographic location. Rhizopus arrhizus is reported to be the most common agent isolated globally, however, other agents such as Apophysomyces is found to be dominant in Asia, Lichtheimia species are found to be dominant in Europe, while Rhizopus homothallicus, Mucor irregularis, and Thamnostylum lucknowense are reported mostly in Asia [32], [33]. The number of epidemiological studies determining the burden of infection is limited, and disease agents have been associated with geographic dispersion. Study findings showed some significant results of mucormycosis coinfection rate reported in certain countries may be due to the frequency of in COVID-19 patients.

From molecular approach, Mucorales bind to the host cell using the endothelium cell receptor glucose-regulated protein 78 (GRP78) [34]–[37]. The pathogenesis of mucormycosis is complicated by the fungus' interaction with the endothelium cells that line the blood vessels. Induced endocytosis allows Rhizopus Oryzae strains to attach to human umbilical vein endothelial cells and invade them. [38]–[40]. Gebremariam et al. found the spore coat protein homolog 3 (CotH3) cell surface protein that binds to GRP78 and examined its role in mucormycosis pathogenesis, where CotH3 acts as a fungal ligand during cell-surface attachment and invasion [41].

The spore coat protein homolog (CotH) cell surface proteins, particularly CotH3, are fungal ligands that enable attachment to the host cell invasion [27]–[29]. Mucormycosis is characterized by vascular infiltration. The fungal ligands for GRP78 were discovered to be Mucorales spore coat protein homologs (CotH). CotH proteins, which were abundant in Mucorales, were absent in noninvasive pathogens [41]. By binding to GRP78, the heterologous development of CotH3 and CotH2 in Saccharomyces cerevisiae was able to gain the ability to penetrate host cells. According to homology modeling research, GRP78 and both CotH3 and CotH2 have structurally compatible interactions [35],[42].

The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) causes the coronavirus disease COVID-19. Understanding SARS-CoV-2 transmission and pathophysiology requires knowledge of the entrance receptor. Early research suggested that the SARS-CoV-2 entrance receptor is angiotensin-converting enzyme 2 (ACE2), the receptor for entrance into the lung epithelial cells [43], [44]–[46]. The cell-surface Glucose Directed Protein 78 (Cs-GRP78) can also act as a multifunctional receptor interacting with many ligands and proteins [47]–[49]. GRP78's substrate-binding domain beta (SBD) has been identified as the binding site for the C480-C488 region within the SARS-CoV-2 spike receptor-binding domain (RBD) [50],[51]. Therefore, CS-GRP78 has been predicted to bind to the SARS-CoV-2 spike protein near the putative host cell receptor ACE2 [52],[53], where GRP78 depleting antibodies blocks the viral entry; hence, GRP78 in association with both ACE2 and the spike proteins would prevent SARS-CoV-2 invasion of the cell [54],[55].

The current work is a computational modeling design of the molecular multi-complex formed by the binding between ACE2 (an entry point in human cells for the coronavirus SARS-CoV-2), RBD (the receptor-binding domain of the coronavirus spike protein), GPR78 (a glucose-regulated protein, which acts as a receptor for host cell invasion of fungi belonging to the Mucorales species, responsible for mucormycosis infection), and CotH3 (a spore coat protein that binds to GPR78 and has a role in the mucormycosis pathogenesis). The assumption is that multiple binding with molecular interactions favors the coinfection between the coronavirus SARS-CoV-2 and the mucormycosis fungi, which would explain why patients with COVID-19 are more likely to develop mucormycosis infection resulting in high morbidity and mortality rates.

2.   Materials and methods

2.1.   ACE2-RBD interactions

The Protein Data Bank (PDB ID: 6M0J) provided the crystal structure of the SARS-CoV-2 spike receptor-binding domain bound to ACE2 (ACE2-RBD). The binding energy of ACE2-RBD (6M0J) was calculated using the end-point free energy calculation approach called the Molecular Mechanics-Generalized Born Surface Area (MM/GBSA) by using HawkDock web server, a tool used for the protein-protein complex structural prediction and analysis [56]. This method is extensively used to estimate the free binding energies and to find the correct binding conformations for protein-protein systems, treating water molecules explicitly is a rigorous technique to account for the solvent effect. The docked structure is then uploaded to PDBePISA for analyzing the interactions of the proteins. Hydrogen bonds, salt bridges, nonbonded contacts, Gibb's free energy of binding, interactive interfaces, tunnels, and pores are all identified in protein complexes. PDBePISA is available at (https://www.ebi.ac.uk/msd-srv/pisa/cgi-bin/piserver?qi=6jpf).

2.2.   GRP78 and ACE2-RBD molecular docking

GRP78/BiP (PDB ID: 5E84) was docked to the spike receptor-binding domain ACE2 of the SARS-CoV-2 crystal structure of the bound with (RBD-ACE2) (PDB ID: 6M0J) using the ClusPro 2 website [57]–[59]. To evaluate the interactions established, the binding energies of the complexes were assessed using the HawkDock webserver's end-point free energy calculation approach MM/GBSA and the PDBePISA tool.

2.3.   CotH3 modeling

The Rhizopus Oryzae CotH3 (RO3G_11882) sequence was obtained from the UniProt database (https://www.uniprot.org). CotH3 was modeled by SWISS-MODEL workspace, a fully automated homology-modeling for the protein structure, using the 5JD9.1 chain A template [60]–[63]. The generated model was then verified using three webservers: PROCHEK, [64],[65], VERIFY 3D [66], and ERRAT [67], available online from the University of California UCLA-DOE LAB using the SAVES v.6.0 webserver (https://saves.mbi.ucla.edu).

2.4.   CotH3 and ACE2-RBD-GRP78 molecular docking

The ClusPro 2.0 webserver for docking protein-protein interactions was used to dock the modeled CotH3 to the RBD-ACE2-GRP78 complex. The binding energies of complexes were calculated using the end-point free energy calculation methodology MM/GBSA of the HawkDock webserver, and the PDBePISA v1.52 tool was used to analyze the interactions formed. The model was examined for SARS-CoV-2 virus variants including the Wildtype PDB ID: 6M0J chain E, Alpha lineage (B. 1.1.7), Beta lineage (B. 1.351), Gamma lineage (P.1), Delta lineage (B.1.617.2), and Delta plus (or Kappa) lineage (B.1.617.2.1) using PyMol software to map and visualize the complexes [68]. The software is available at (https://pymol.org).

3.   Results

The binding energy of ACE2-RBD (6M0J) was first calculated using the MM/GBSA of the HawkDock webserver. The PDBePISA tool was then used to evaluate the expected interface features and key atoms between the docked structures, such as hydrogen bonds and salt bridges, using the PDBePISA tool. The predicted salt bridges and hydrogen bonds are listed in Table 1. The number of hydrogen bonds and salt bridges is used to assess the likely stability of the interface.

ClusPro 2.0 was then used to dock GRP78 to ACE2-RBD. The resulting complexes were ranked by the binding energies calculated using MM/GBSA where the complex that had the highest binding energy was selected. The results are displayed in Figure 1. The binding energies between GRP78 and ACE2 with RBD within the ACE2-RBD-GRP78 complex were calculated using MM/GBSA. Table 1 Table 1 displays the results of the PDBePISA for the salt-bridge and hydrogen bond interface interactions.

The CotH3 model of Rhizopus Oryzae was created by the SWISS-MODEL workspace using the 5JD9.1 chain A template (18.68% sequence similarity). The model is shown in Figure 2A. According to the Ramachandran plot, 96.6% of the model residues are within the allowed regions, as shown in Figure 2B. The VERIFY tool assessed approximately 81.5% of the residues to have a 3D-1D score ≥ 0.2, and the overall ERRAT quality factor was 86.12%.

ClusPro 2.0 was used again to dock the CotH3 model to the ACE2-RBD-GRP78 complex. The binding energies of CotH3 and ACE2, RBD, and GRP78 were calculated using MM/GBSA. The complex that had the highest binding energy was selected. Results are ranked by the binding energies and shown in Figure 3. PDBePISA was then used to predict the number of the salt bridges and hydrogen bonds of the interacting interface. Results are listed in Table 2 and the interacting sites are displayed in Figure 4. Approximately 2961 mutations in Spike protein have been found based on the genome sequences of the SARS-CoV-2. We obtained information on all S protein mutations found in the RBD domain. Table 3 displays the binding energies of CotH3 and SARS-CoV-2 RBD for the lineages examined including Alpha, Beta, Gamma, Delta, and Delta plus variants.

4.   Discussion

The findings of GRP78 binding to the RBD are consistent with Ibrahim et al. study findings [50], which presented that GRP78 binds to 4 regions with the Spike protein of the SARS-CoV-2 virus, among which regions III (C391-C525) and IV (C480C488) show stronger Spike RBD affinity [46]. The results show that the substrate-binding domain β of GRP78 weakly binds to region IV of WT RBD (C480-C488) (−18.5 kcal/mol) forming three hydrogen bonds with PRO479, CYS480, and LYS481, as shown in Table 1. Figure 1 and Figure 3 show that GRP78 and ACE2 bind to the spike RBD surface in roughly the same sites and bind to each other. This is consistent with the findings of Aguiar et al. [69] who showed that the GRP78-binding location overlaps with the ACE2-binding location, albeit the residues engaged in the interactions may differ somewhat. The binding energy between the RBD and ACE2 is increased by around 45%, from −60.55 kcal/mol to −87.5 kcal/mol due to GRP78. In the ACE2-RBD complex, RBD forms 11 hydrogen bonds and two salt bridges with ACE2, whereas in the GRP78-RBD-ACE2 complex, the complex forms 14 hydrogen bonds and four salt bridges with ACE2, as shown in Table 1. Therefore, GRP78 stabilizes the binding of RBD and ACE2, which improves the successful entry of the virus.

The results of CotH3 binding to the ACE2-RBD-GRP78 complex show that while CotH3 weakly attaches itself to GRP78 (−2.63 kcal/mol) and repels from ACE2 (+21.15 kcal/mol), it strongly binds to RBD (−116.91 kcal/mol) as it forms 16 hydrogen bonds and doubles the number of salt bridges with the spike RBD, as listed in Table 2. CotH3 shows a similar affinity for the SARS-CoV-2 virus variants of RBD with some preference for the Delta variants as shown in Table 3. This means that the fungus can enter the cell via the spore coat protein attached to the SARS-CoV-2 RBD spike, which would explain the high coinfection rate of mucormycosis in COVID-19 patients.

5.   Conclusions

Mucormycosis, known as black fungus, is a life-threatening fungal infection. The growing number of case mortality and morbidity reported globally during the second wave of the Corona Virus Disease 2019 infection (COVID-19) suggested a high coinfection rate of mucormycosis among COVID-19 patients, where the infection rates of Mucormycosis in COVID-19 patients have been increased.

Current work presents a hypothesis of the molecular interactions that favor the coinfection between the coronavirus SARS-CoV-2 and the mucormycosis fungi. The present research work explains the observed increased morbidity and mortality coinfection in COVID-19 patients with mucormycosis by using the proposed computational molecular model of the fungus spore coat protein that binds to the spike protein of the SARS-CoV-2 virus. The SARS-CoV-2 RBD complex with ACE2 and GRP78 binds to the CotH3 ligand of mucormycosis, carrying it into the host cells during COVID-19 infection. The findings from this study explain that the COVID-19 infection generates an appropriate environment for the spread of Mucorales, resulting in a high infection rate of mucormycosis. The outcome from this study suggests that COVID-19 patients with severe symptoms be evaluated for mucormycosis infection when comorbidities are present, and further research is needed to evaluate the potential association between both infections.