Screening of potential antigens from whole proteome and development of multi-epitope vaccine against Rhizopus delemar using immunoinformatics approaches

Abstract Mucormycosis is a deadly fungal disease mainly caused by Rhizopus oryzae (strain 99–880), also known as Rhizopus delemar. Previously, mucormycosis occurs in immunocompromised patients of diabetes mellitus, cancer, organ transplant, etc. But there was a drastic increase in mucormycosis cases in the ongoing COVID-19 pandemic. Despite several available therapies and antifungal treatments, the mortality rate of mucormycosis is about more than 50%. Currently, there is no vaccine available in the market for mucormycosis that urgently needs to develop a potential vaccine against mucormycosis with high efficacy. In the present study, we have screened 4 genome-derived predicted antigens (GDPA) through sequential filtration of the whole proteome of R. delemar using different benchmarked bioinformatics tools. These 4 GDPA along with 4 randomly selected experimentally reported antigens (ERA) were sourced for prediction of B- and T- cell epitopes and utilized in designing of two potential multi-epitope vaccine candidates which can induce both innate and adaptive immunity against R. delemar. Besides these, comparative immune simulation studies and in silico cloning were performed using L. lactis as an expression system for their possible uses as oral vaccines. This is the first multi-epitope vaccine designed against R. delemar through systematic pipelined reverse vaccinology and immunoinformatic approaches. Although the wet-lab based experimental validation of designed vaccines is required before testing in the preclinical model, the current study will significantly help in reducing the cost of experimentation as well as improving the efficacy of vaccine therapy against mucormycosis and other pathogenic diseases. Communicated by Ramaswamy H. Sarma


Introduction
Mucormycosis, also known as black fungus, is an invasive fungal disease caused by fungi belonging to the order Mucorales. The term mucormycosis was first proposed by Baker in 1957(Baker, 1957. About 11 genus and approximately 27 species of Mucorales are responsible for mucormycosis in humans. Out of these, mucormycosis is primarily caused by Rhizopus oryzae in humans which covers nearly 70% of total mucormycosis cases (Jeong et al., 2019;Prakash & Chakrabarti, 2019). The disease mucormycosis predominantly occurs in immunocompromised hosts because of diabetes mellitus, iron overload, neutropenia, major trauma, burns, treatment with corticosteroids, bone marrow transplant, organ transplant, hematopoietic stem cell transplantation (HSCT), malignant hematologic disorders, deferoxamine therapy in patients receiving hemodialysis, COVID-19 or severe immunosuppression but can also develop in immunocompetent hosts. It may affect adults as well as children and its mortality rate is about 50% (Castillo et al., 2018;Ibrahim et al., 2003;2012;Ma et al., 2009;Spellberg et al., 2005;Sugar, 2005). Different forms of mucormycosis are rhino-cerebral mucormycosis, pulmonary mucormycosis, gastrointestinal mucormycosis, renal mucormycosis, disseminated mucormycosis, and cutaneous mucormycosis which mainly occur due to infection of sporangiospores in nasal and oral during inhalation, ingestion of contaminated food, infection in conjunctiva mucosa, and traumatic inoculation. It cannot be transferred from human to human (Paduraru et al., 2016;Prakash & Chakrabarti, 2019;Wilson et al., 2019).
The 45 Mbp genome sequence of a clinical Rhizopus oryzae strain 99-880 (recently reclassified as Rhizopus delemar) was first published in 2009 which comprises about 14,000 predicted protein-coding genes and a high proportion of repetitive transposable elements. It is thought to have undergone a whole-genome duplication (WGD) event followed by massive gene loss (Ma et al., 2009). The WGD resulted in the expansion of gene families related to fungal-specific cell wall synthesis, pathogen virulence, and signal transduction. This WGD of fungi might be responsible for the rapid growth and tremendous adaptive capacity in the host, adaptation to the environment, immune escape mechanism, antifungal resistance, etc. (Kontoyiannis & Lewis, 2015;Ma et al., 2009). During WGD of R. oryzae, there was an expansion of 23 chitin synthases (CHS) and 24 chitin deacetylases (CDA) gene families that might be involved in the synthesis of chitin and chitosan, respectively. Chitosan or chitin are the main components of the fungal cell wall which is involved in immune escape mechanism and can be used as a target for antifungal development, diagnosis, and therapeutics (Ma et al., 2009;Sephton-Clark et al., 2018).
Usually, when R. delemar spores enter inside the host, it must be destroyed by phagocytosis. But in the case of immunocompromised conditions or impaired phagocytic function, these spores actively change into hyphae and can cause mucormycosis. Scientific literature has revealed that overexpression of glucose-regulated protein 78 (GRP78) receptors on the cell surface is associated with endoplasmic reticulum stress conditions (diabetic/immune suppression). During the immunocompromised condition, the CotH (mostly CotH3) surface protein of R. delemar interacts with GRP78 receptors of the host and invades endothelial cells (Gebremariam et al., 2019). After crossing of endothelial cells, the resting spores of R. delemar adhere to type IV collagen and laminin protein of the basement membrane. These resting spores of R. delemar convert into swollen spores followed by hyphae. These hyphae further invade the vasculature and enter the lumen of blood vessels where they form a thrombus and lead to ischemia, sustained hypoxia, and necrosis (Ghuman & Voelz, 2017;Liu et al., 2010).
Early clearance of R. delemar infection is associated with phagocytes and innate immune response. The surface of fungal spores and hyphae expresses pathogen-associated molecular patterns (PAMPs) which are recognized by patternrecognition receptors of phagocytes such as toll-like receptors (TLR) and induce an innate immune response (Romani, 2011). The TLR2 can recognize the PAMPs of fungal hyphae to induce innate immune system and inflammatory responses such as stimulation of macrophage to release TNF-a and cytokines which eventually help in the clearance of fungi (Chamilos et al., 2008;Figueiredo et al., 2011;Park & Voigt, 2014;Roilides et al., 2012). Both platelets and monocytes expressed TLR2 and TLR4 on their surface which might be associated with fungal defense (Bruserud, 2013). Natural killer (NK), neutrophils, and dendritic cells (DC) are also associated with innate immune response activation to provide protection against mucormycosis. However, the adaptive immune response is induced by activated DC (Ghuman & Voelz, 2017). Besides these, helper T (Th) cells mediated inductions of a cytokine such as IL-4, IL-10, and interferon (IFN) are helpful in the clearance of R. delemar (Potenza et al., 2011).
Before the COVID-19 pandemic, the rate of mucormycosis varies from 0.005-1.7 cases per 1,000,000 population globally. But in India, it is approximately 0.14 cases per 1000 population, which is about 80 times greater than that of the worldwide cases (Chander et al., 2018;Maini et al., 2021). In the recent past, mucormycosis infections were abruptly increased worldwide due to the COVID-19 pandemic (John et al., 2021;Johnson et al., 2021;Veisi et al., 2021;Waizel-Haiat et al., 2021). In India, according to the official data published in a leading newspaper, large number of cases/deaths from Maharashtra (1500/90), Gujarat (1163/61) , Madhya Pradesh (575/31), Haryana (268/8), Delhi (203/1), Uttar Pradesh (169/ 8), Bihar (103/2), Chhattisgarh (101/1), Karnataka (97/0) and Telangana (90/10) have been reported for mucormycosis in people who have suffered from COVID-19 till 20 th May 2021 (Hindustan Times, 2021, Published on May 21, 2021, 07:09 AM IST). Different research groups all over the world reported that severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infections enhance the risk of mucormycosis in patients (Maini et al., 2021;Mehta & Pandey, 2020;Revannavar et al., 2021). Despite several therapies such as antifungal therapy (monotherapy or combination therapy), reversal of immunosuppression (an increase of neutrophils and discontinuation of chemotherapy), and surgical resection of necrotic tissues, etc., are helpful for the management of mucormycosis but it is still a life-threatening disease (Sipsas et al., 2018) which might be prevented through vaccination strategy. Currently, there is no licensed vaccine available in the market against mucormycosis , WHO Emergency Report, 2021. But, because of the high morbidity and mortality of mucormycosis, there is an urgent need to develop an efficient vaccine that can induce a broad immune response to prevent infections in highly immunocompromised as well as immunocompetent hosts (Spellberg, 2007).
The present study utilized the reverse vaccinology and immunoinformatics approaches for designing of new vaccine candidates against mucormycosis involving whole protein sequences of R. delemar. The study design involved the screening of genome-derived predicted antigens (GDPA) and randomly selected experimentally reported antigens (ERA) for the prediction of B-and T-cell epitopes, designing of vaccine constructs, and evaluation of their immunological properties through docking and immune simulation.

Retrieval of protein sequences of R. delemar
The entire protein sequences of R. delemar were retrieved from the FungiDB database (Stajich et al., 2012) for the screening of GDPA. However, the ERA was randomly selected through a comprehensive review of the literature and their protein sequences were retrieved from the NCBI database (NCBI Resource Coordinators, 2018).

Prediction of T cell epitopes
The IEDB based MHC class I and MHC class II binding tools were used for the prediction of MHC class I and MHC class II binding epitopes for each of the 4 GDPA and 4 ERA with the selection of all the alleles and consensus method at respective thresholds IC 50 50 nM and percentile rank 1 (Fleri et al., 2017). Since most of the MHC class I molecules prefer to bind 9-mer peptides due to the closed binding cleft therefore 9-mer length of the peptide was selected for MHC class I epitope prediction (Nielsen & Andreatta, 2016) while in the case of MHC class II molecule the default value of peptide length 15-mer was used (Wang et al., 2008;. MHC class II proteins usually accommodate peptides of variable length ranging from 3-25 residues (mostly 15-mer) in their open binding groove with 9-mer core motifs (Nielsen & Andreatta, 2016;Wang et al., 2008;. Furthermore, the immunogenicity of epitopes was predicted using IEDB based MHC class I Immunogenicity tool (Calis et al., 2013). The predicted population coverage (PPC) value of epitopes was predicted through IEDB based population coverage analysis tool (http://tools.iedb.org/ population/) against the worldwide population. Based on the criteria of highest PPC value, VaxiJen score, and immunogenicity score, a minimal T cell epitope ensemble with maximum PPC value was selected for both GDPA and ERA excluding the epitopes with the same set of HLA binding alleles as detailed in a previous study by Pritam et al. (2019). The toxicity of epitopes was also predicted using the ToxinPred tool .

Designing of vaccine candidates
A simple protocol reported in our previous study of Pritam et al. (2020) was used to design two new vaccine candidates namely ERAMV1 and GDPAMV1 corresponding to the identified epitope ensemble of ERA and GDPA using suitable linkers EAAAK, GGGS, and GPGPG (Chen et al., 2013;Dong et al., 2020;Meza et al., 2017). The N-terminal sequence of the vaccine candidate was initiated with an EAAAK linker followed by an adjuvant, MHC class I epitopes, MHC class II epitopes, and linear B cell epitopes. The linker GGGS was used for connecting MHC class I epitopes while the GPGPG linker was used to join MHC class II epitopes. Besides these, GPGPG linker was also used for joining of MHC class I epitopes and MHC class II epitopes together as well as linking of linear B-cell epitopes. The Cholera toxin B subunit (A: UniProt accession no. AIE88420.1) was used as an adjuvant against TLR2 (PDB ID: 2Z7X) (Jin et al., 2007;Pritam et al., 2020).

Prediction of discontinuous B-cell epitopes
The discontinuous B-cell epitopes from the 3 D structure of both designed vaccine candidates ERAMV1 and GDPAMV1 were predicted using IEDB based Ellipro tool at threshold 0.6 (Ponomarenko et al., 2008).

MHC class I processing analysis of designed vaccine
The MHC class I processing of designed vaccines was performed by using the MHC class I processing predictions tool of IEDB at default parameters (http://tools.iedb.org/processing/). All the predicted MHC class I processing epitopes for MHC class I was screened at threshold IC 50 500 nM. For the analysis of MHC class I processing of epitopes, we have used a threshold of 1, 0.3, and 1.5 for proteasome score, TAP score, and processing score, correspondingly (Tenzer et al., 2005).
Molecular docking of designed vaccine candidates with TLR2 and calculation free energy using MM/ GBSA approach The molecular docking and interaction analysis of TLR2 with designed vaccine candidates and co-crystallized ligands (positive controls) were performed using Cluspro 2.0 (https:// cluspro.bu.edu/home.php) (Kozakov et al., 2017) and discovery studio visualization tools, respectively. Three positive control ligands viz. Escherichia coli heat-labile enterotoxin type IIB B-pentamer (C1; PDB ID: 1QB5) and Pam3CSK4 (C1A; PDB ID: 2Z7X) and SSL3 (C1B; PDB ID: 5D3I) were used in the docking studies for comparative evaluation of binding affinity with ERAMV1 and GDPAMV1 (Koymans et al., 2015;Zhong et al., 2015). The free energy of interaction between the TLR2 and designed vaccine candidates (ERAMV1 and GDPAMV1) was estimated by using the molecular mechanics generalized-born surface area (MM/GBSA) approach utilizing HawkDock server (Weng et al., 2019). In addition, the HADDOCK (https://wenmr.science.uu.nl/haddock2.4/submit/1) webserver was also used for molecular docking at default parameters (van Zundert et al., 2016). After analysis of interaction sites of ligand and receptor of 2Z7X and 5D3I complex, we have selected a region 291 to 355 position of TLR 2 residues as highly probable binding sites.

Molecular dynamics analysis of TLR2-designed vaccine complex
The molecular dynamics study of docked models of TLR2 and designed vaccine candidates were performed using a normal mode analysis (NMA) approach with the iMODS web server (http://imods.chaconlab.org/). The NMA predicts the collective motion of atoms of a protein complex (L opez- Blanco et al., 2014).

Disulfide engineering of designed vaccine candidates
The tool Disulfide by Design 2.0 was used for disulfide engineering of designed vaccine candidates. It predicts the amino acid residue that can be replaced with cysteine residues for the formation of a disulfide bond.
In silico cloning of designed vaccine candidates along with the prediction of cDNA, codon optimization, and prediction of mRNA secondary structure The JCAT server was used for the prediction of cDNA and codon optimization of designed vaccines (Grote et al., 2005). For in silico cloning of cDNA of designed vaccines, the SnapGene software was used, and obtained cDNA was inserted in plasmid vector pIL1 (Gene bank accession number: HM021326) at FspI (6006) restriction site (Hayes et al., 1990). The model organism Lactococcus lactis (strain IL1403) was used as an expression host and avoided the rho-independent transcription terminators, prokaryotic ribosome binding sites, and cleavage sites of restriction enzymes. The optimized cDNA was converted into mRNA by using a transcription and translation tool at Cybertory (http://biomodel. uah.es/en/lab/cybertory/analysis/trans.htm). Finally, the RNA Folding Form version 2.3 of Mfold tool available at the UNAFold web server (http://www.unafold.org/mfold/applications/rna-folding-form-v2.php) was used for prediction of the secondary structure of mRNA at default parameters including folding temperature at 37 C (Zuker, 2003).

Whole proteome-based screening of antigens
A total 17459 protein sequences of R. delemar were found in the FungiDB. The sequential filtration from 17459 protein sequences through SignalP, PredGPI, VaxiJen2.0, and TMHMM tools lead to the screening of 684, 64, 43, and 42 proteins, correspondingly (Supplementary Table 1). Out of screened 42 GPI-anchored antigenic proteins with signal peptides and 1 transmembrane helix (TMH) region, we have selected the top 4 GDPA viz. RO3G_04553, RO3G_02456, RO3G_14327, and RO3G_05587 for further study purely based on VaxiJen2.0 score. Certainly, a higher antigenic score of a protein depicts its high probability for possible involvement in the activation of the immune system. Proteins with 1 TMH region were selected because a transmembrane protein containing >1 TMH region can cause difficulty in the process of cloning, expressing, and purification (He et al., 2010). Signal peptides are also the peculiar signals for the translocation of secretory proteins across membranes in both the eukaryotic and prokaryotic systems. Precisely, in the eukaryotic system, several proteins have been found to link with the extracellular leaflet of the plasma membrane, which transports a GPI-anchor associated with the C-terminal residue later during the process of proteolytic cleavage occurring at omega-site (Guy et al., 2018;Singh et al., 2020). These surface-associated-GPI anchored proteins are found to be associated in process of immune evasion and cyto-adhesion; hence these GPI-anchored proteins could be a suitable target to be considered as vaccine candidates. Therefore, in the present study, these 4 GDPA were used as potential candidates for the screening of B-and T-cell epitopes to be used in designing a multi-epitopes vaccine against mucormycosis. A similar vaccine based on the surface protein of Candida albicans elicited high levels of antibodies (IgG, IgM) and showed a high degree of protection in mice (Xin et al., 2012). Hence, to compare the predictability of present bioinformatics approaches for screening of genome-derived vaccine candidates, a literature survey based set of 4 randomly selected ERA viz. Hypothetical protein RO3G_11882 (accession no. EIE87171.1), hypothetical protein RO3G_17119 (accession no. EIE92521.1), hypothetical protein G6F43_000212 (accession no. KAG1057968.1), and endopeptidase (accession no. EIE89256.1) were obtained from the NCBI gene bank database (Table 1) (Carroll et al., 2017;Gebremariam et al., 2019;Sassone-Corsi et al., 2016;Sato et al., 2017;Sircar et al., 2012). All the antigens belong to Rhizopus delemar strain RA 99-880 except for hypothetical protein G6F43_000212 (Rhizopus delemar strain GL54). Amongst these, the hypothetical protein RO3G_11882 is also known as CotH3 which is directly associated with host cell interaction and invasion. Additionally, the CotH (CotH3) protein is highly targeted for therapeutic purposes, such as in the development of an anti-CotH antibody (Gebremariam et al., 2019).

Prediction of linear B-cell epitopes
The amino acid residues exposed on the surface of antigenic proteins that can interact with the antibodies are known as B-cell epitopes. A total of 28 and 24 linear B-cell epitopes were predicted from a set of 4 GDPA and 4 ERA using the BCPREDS tool, respectively (Supplementary Table 2). Out of these, only the top 2 antigenic B-cell epitopes from each protein were utilized in multi-epitopes vaccine designing ( Table 1).

Prediction of T cell epitopes
T-cells are classified into helper T cells (Th) and cytotoxic T cells (Tc) which contain CD4þ and CD8þ membrane protein, correspondingly. The respective co-receptors CD4þ and CD8þ interact with the peptide-major histocompatibility complex (p-MHC) II and I ligands along with T-cell receptors (TCRs) and facilitate cell-mediated immunity against fungi (Artyomov et al., 2010). The murine model studies provided evidence for Th and Tc associated control against fungal disease through induction of cytokines (IFN-c, IL-10, IL-17) and cytotoxic activity, correspondingly (Antunes et al., 2018;De Luca et al., 2012). Keeping these facts in mind, the present study predicted a total of 52 and 47 MHC class I and class II binding epitopes from 4 ERA as immunogenic epitopes, respectively. While in the case of GDPA, a total of 30 immunogenic epitopes were predicted against each MHC class I and MHC class II molecules (Supplementary Table 3). The number of above-mentioned epitope sets was further reduced to ensemble epitope set of 16 and 7 for ERA and 12 and 4 for GDPA based on the highest PPC value, VaxiJen score, immunogenicity score as the strategy adapted in our previous study of Pritam et al. (2020) (Table 2). The respective epitope ensemble PPC against the worldwide population of GDPA 96.16%, 79.29%, and 99.20% were found significantly higher than the ERA 83.78%, 51.73%, and 92.17% for individual MHC class I, MHC class II and combined (both MHC class I and II) (Figure 2). These ensemble epitopes were also found non-toxic and hence used in designing multi-epitope vaccine candidates.

Prediction of cytokines responses of ensemble epitope
Induction of adaptive immune response mediated by T-cells can produce different cytokines such as interleukin-2 (IL-2), IL-4, IL-5, IL-7, IL-10, IL-13, IL-1b, interferon-c (IFN-c), and tumor necrosis factor (TNFa) which are found to be associated with clearance of R. oryzae (Dos Santos et al., 2020;Page et al., 2018;Potenza et al., 2011). These cytokines can be used as immunotherapy against mucormycosis because the treatments of IL-2/IL-7 elicit IL-13, IL-5, TNFa, and IL-10 responses. Also, IL-2 is associated with the stimulation of Tcells, and IFN-c is associated with hyphal damage of fungi and phagocytes activation that can protect R. oryzae (Castillo et al., 2018; Potenza et al., 2011). These facts provided sufficient foundation for in silico evaluation of ensemble epitope corresponding to GDPA and ERA for cytokines (IL-4, IL-10, and IFNc) response induction (Supplementary Table 4). Fortunately, the epitopes E18, and E23 of ERA were predicted as IL-4, IL-10, and IFN-c inducer while E3 and E8 were predicted only as IL-4 and IFN-c inducer. However, the epitope G6 of GDPA was also predicted as IL-4 and IFN-c inducer. Since the IL-17 is also helpful in the clearance of R. delemar, the present study developed a machine learning-based prediction model for the evaluation of IL-17 induction response prediction (Potenza et al., 2011). Initially, for the development of the prediction model, 1255 positive and 6803 negative epitopes associated with T-cell assay for IL-17 induction were retrieved from the IEDB database. The 9-mer binding core motifs associated with IL-17 were identified through sequence alignment using GibbsCluster 2.0 . The identified 9-mers binding core motifs of positive (744) and negative (3042) epitopes dataset were used for training and evaluation of the IL-17 inducing prediction model. In this regard, a random method was used to create training and validation subsets along with blosum amino acid numerical encoding. The training of ANN method using NNAlign-2.0 tool involved the following parameters viz. motif length (9), number of hidden neurons (10), number of training cycles (500), number of random seeds (4), number of networks in the final ensemble (20), and number of folds for cross-validation (5). The prediction performance of the developed ANN model was evaluated using three different metrics viz. rootmean-squared error (RMSE), Pearson correlation coefficient (PCC), and Spearman rank coefficient (SRC). The RMSE reflects the error between the mean of the experimental and predicted  (Mukaka, 2012;Veerasamy et al., 2011). Finally, the own developed model predicted 14 and 9 epitopes of ERA and GDPA as IL-17 inducers at threshold > 0.5 (Supplementary Table 4).

Designing of multi-epitope vaccine candidates
In the last few years, immunoinformatics approaches are being very frequently used for the design and development of multi-epitope based vaccines against various diseases such as COVID-19, Leishmaniasis, Schistosomiasis, Malaria, Diarrhea, Influenza A, Malignant melanoma, Dengue, Fascioliasis, etc. because of convenience, low cost, reduced labor and/or risk of failure of the vaccine in the wet laboratory (Ali et al., 2017;Behbahani et al., 2021;Kaba et al., 2018;Kazi et al., 2018;Rahmani et al., 2019;Singh et al., 2020;Yazdani et al., 2020;Zeb et al., 2021). In this study, we reported the design of two new vaccine candidates namely ERAMV1 and GDPAMV1 utilizing the respective ensemble epitope of ERA and GDPA along with adjuvant and suitable linkers EAAAK (L1), GGGS (L2), and GPGPG (L3) ( Table 3). A linker used in multi-epitope vaccine design provides a stable structure, facilitates the antigen-presenting process, enhances immunogenicity, and prevents the arrangement of junctional epitopes (Livingston et al., 2002;Nezafat et al., 2016;Yang et al., 2015). The linker L1 was used to provide inflexibility which averts the assemblage of adjuvant with other domains of the vaccine while linker L2 was used to increase bendability and immune response. The universal linker L3 was used to facilitate the antigen processing with a proteasome that ultimately enhances immunogenicity Yang et al., 2015;Yano et al., 2005). To further enhance the overall immunogenicity (both innate as well as adaptive) and degree of protection against fungi, the cholera toxin B (CTB) adjuvant was used as a TLR2 activating agent. The TLR2 plays an important role in the activation of autophagy in macrophages which can protect against fungi during the early stage of infection (Wu et al., 2019). The interaction of TLR2 and fungal hyphae can induce the TLR2 mediated immune response by activating innate immune cells such as neutrophils. These activated innate immune cells (neutrophils) can further cause damage to the hyphae of fungus (Ghuman & Voelz, 2017). In the case of R. oryzae, TLR2 also plays a significant role to prevent fungal infection by inducing a pro-inflammatory response (Chamilos et al., 2008). The CTB adjuvant is cost-effective and safe, as well as facilitates the immune induction mediated by Th with a high level of antibody titration. The CTB has been evaluated as an adjuvant in several scientific studies for the development of an effective vaccine against the influenza virus and showed significant induction of both cellular as well as humoral immune responses (Lei et al., 2015;Stratmann, 2015).

Analysis of designed vaccine candidate properties based on amino acid sequence
A potential vaccine should have the properties like hydrophilic in nature, stability and durability (Bastola & Lee, 2019;Guruprasad et al., 1990). In the past, huge numbers of hydrophilic vaccines were reported by different research groups that induce innate, humoral, and cell-mediated immune responses (Ali et al., 2017;Solanki et al., 2019). However, in the last few decades, several hydrophobic protein-based vaccines have been also evaluated for their safety and found to induce immune responses including antibody response (IgG and IgA) and CD4þ T cell-dependent IFN-c response (Azuar et al., 2021;Langley et al., 2018;Schepens et al., 2015;Torrey et al., 2020). The hydropathy analysis (based on the GRAVY score) of designed vaccines showed that ERAMV1 is highly negative while GDPAMV1 was slightly positive which reflects their hydrophilic nature and slightly hydrophobic nature, correspondingly. The instability index value of both ERAMV1 and GDPAMV1 were found < 40 and showed 1 hour in vivo half-life hence implicating enough stability (Table 4). The obtained results of tools predicting antigenicity such as VaxiJen 2.0, ANTIGENpro, Protein antigenicity prediction, and Secret-AAR indicate that both the designed vaccines were antigenic in nature (Table 4). The protein function of designed vaccines predicted by DeepGOPlus validated their potential role in immune induction through the property of killing of cells of another organism instead of the host organisms. Both the designed vaccines were found soluble as predicted by Protein-Sol, CamSol, and SOLPro, and secondary structure analysis (alpha helices, b-strands, and coils) showed the presence of regular structures like alpha-helices, b-strands, and coils for proper folding, and stability (Berg et al., 2002;Buchan & Jones, 2019;Ding et al., 2014;Lovell et al., 2003;Rost, 2001) (Figure 3; Table 4). Furthermore, no significant similarity was found against Homo sapiens for both designed vaccines.

Prediction of the tertiary structure of designed vaccine candidates
The tertiary structure of the designed vaccines ERAMV1 as well as GDPAMV1 showed the proper three-dimensional (3 D) folding ( Figure 4). Ramachandran plot was frequently used by many bioinformaticians for the analysis of the accuracy of the tertiary structure of the protein (DasGupta et al., 2015). Thus, the quality of protein vaccine 3 D structures was evaluated using Ramachandran plot through Ramachandran plot web server. For a good quality model of a protein, over 88% of amino acid residues are expected to be found in the most favored regions. The amino acid residues of both the designed vaccines ERAMV1 (89.9%) and GDPAMV1 (92%) were found in most favored regions which reflects the good quality of stable protein 3 D structures ( Figure 5). Further, the result obtained by Verify3D showed 90.08% and 84.16% residues of ERAMV1 and GDPAMV1 have average 3 D-1D score > ¼ 0.2, respectively, and passed the tertiary structure verification. The RMSD, and TM scores of refined models of ERAMV1 (2.482 and 0.935) GDPAMV1 (1.531 and 0.967) indicated that 3 D structures were properly folded.

Prediction of discontinuous B-cell epitopes
A total of 4 and 2 sets of discontinuous B-cell epitopes were predicted from the 3 D structure of designed vaccines ERAMV1 and GDPAMV1, respectively at a threshold above 0.6 (Supplementary Table 5; Figure 6).

MHC class I processing analysis of designed vaccine
Whenever an antigenic protein enters into the antigen-presenting cell it undergoes proteasomal degradation to produce small peptides. These peptides are further transported by TAP (transporter associated with antigen processing) to endoplasmic reticulum lumen where they can be uploaded to MHC class I molecule and ultimately expressed on the cell surface for stimulation of CD8þ cytolytic T cells (CTL) response (Vigneron et al., 2017). After analysis of MHC class I processing of ERAMV1 and GDPAMV1, we have found that most of the MHC class I epitopes of designed vaccines were obtained as CTL epitopes after proteasomal degradation. In the case of ERAMV1, epitopes E2, E4, E5, E6, E7, E13, E14, E16, E19, E21, and E22 were predicted as CTL epitopes. While in the case of GDPAMV1, epitopes G2, G5, G7, G10, G11, and G14 were predicted as CTL epitopes (Supplementary  Table 6).

Molecular docking of designed vaccines with TLR2 and calculation of MM/GBSA based binding free energy
The molecular interaction studies between the designed vaccines along with positive control (C1) and TLR2 were evaluated through the protein-protein docking protocol of the ClusPro 2.0 tool. Both the designed vaccines showed lower docking energy compared to the positive control (C1) that indicates good molecular interaction between designed vaccines and TLR2 and that generates the possibility to induce TLR2 associated immune response (Table 5; Figure 7). A set of common interacting residues of TLR2 were also found between designed vaccines and controls (C1, CA1, and CB1) ( Table 6). The SER-2 residue of control C1 interacts with PHE-325 and ASP-327 residues of TLR 2 of complex crystal structure (PDB ID: 2Z7X). Further, the estimated binding free energy by using MM/GBSA approach utilizing HawkDock server of molecular complexes of TLR2-ERAMV1 (-14.77 Kcal/ mol) and TLR2-GDPAMV1 (-4.54 Kcal/mol) revealed a strong interaction of TLR2 and designed vaccines ERAMV1 and GDPAMV1. In addition, the Z-Score of HADDOCK 2.4 also indicates the good interaction between the designed vaccines and TLR2 compared to control (C1). The more negative Z-score indicates the better binding affinity (Rostaminia et al., 2021;van Zundert et al., 2016). Thus, the binding affinity of GDPAMV1 with TLR2 was found to be better than ERAMV1 and control (C1) ( Table 5).

Molecular dynamics analysis of the docked complexes of TLR2 and designed vaccines
Molecular dynamics study is crucial for evaluating the stability of the protein-protein complex, which can be determined by comparing necessary protein dynamics to their normal modes (W€ uthrich et al., 1980). The stability of the proteinprotein complex can be predicted by molecular dynamics analysis associated with normal mode analysis (NMA). Thus, the mobility and stability studies of protein-protein complexes of the designed vaccines and TLR 2 were evaluated through NMA (Figure 8). The 3 D interaction model of EAMV1 and GDPAV1 revealed the direction of each amino acid residue by arrows and the length of the arrow indicates the corresponding degree of mobility (Figure 8a and b). The resulting output also provided the deformability (Figure 8c and d), mobility (Figure 8e and f), eigenvalue (Figure 8g and h), variance map (Figure 8i and j), covariance matrix ( Figure  8k and l), and elastic network (Figure 8m and n). The relative amplitude of the atomic displacements around the equilibrium confirmation was indicated by the value of NMA-B-factors (mobility). The gradient of the atomic displacements at every atomic position summed over all modes was presented by deformability. Higher value reflects the flexibility of complex structure while lower value indicates rigidity. The obtained higher and lower values of maximum mobility and deformability for ERAMV1 (9.972E þ 00, 2.579E-06) indicated more flexible regions compared to GDPAMV1 (4.381E þ 00, 4.196E-06). The eigenvalue indicates the motion of stiffness associated with each normal mode. For easy deformation, the eigenvalue should be lower which reflects that lower energy is required for better results with deformation of complex structure. The respective eigenvalues 1.68e-06 and 6.58e-06 of ERAMV1 and GDPAMV1 reflect the stability of the TLR2-ERAMV1/PAMV1 complexes. The individual and cumulative variances associated with each normal mode were inversely related to the eigenvalue. The covariance matrix presented the association of pairs of residues. Correlated, uncorrelated, and anti-correlated relation was presented in red, white, and blue color, correspondingly. In an elastic network graph, each dot represents one spring characterized by pairs of atoms (Azim et al., 2019).

Disulfide engineering of designed vaccines
Disulfide engineering is an important component of multiepitopebased vaccine designing to enhance functional characteristics, stability and assists molecular dynamics of protein (Craig & Dombkowski, 2013). It works by increasing the free energy of the denatured state and reducing conformational entropy which ultimately increases the stability of folded protein conformation (Flory, 1956 Table 7). These predicted disulfide bonds can significantly enhance the stability of the designed vaccines (Mugunthan & Harish, 2021). In silico cloning of designed vaccines and prediction of mRNA secondary structure The results from the JCAT web server provide the prediction and codon optimization of cDNA of designed vaccines ERAMV1(1965 bp) and GDPAMV1 (1629 bp), correspondingly. For the in silico cloning of cDNA of designed vaccines, the obtained cDNA was inserted into the plasmid vector pIL1 (Gene bank accession number: HM021326) at FspI (6006) restriction site using SnapGene software (Grote et al., 2005;Hayes et al., 1990). The obtained reliable values of respective codon adaptation index (CAI) and GC-Content for ERAMV1 (0.94: 55.77%) and GDPAMV1 (0.97: 55.06%) were found above the optimal range i.e., 0.9-1.0 and 30%-70%, correspondingly (Grote et al., 2005; Urrutia-Baca et al., 2019). Accordingly, both the designed vaccines can be easily expressed in the L. lactis expression host. After insertion of cDNA into a plasmid, the size of obtained recombinant DNA of ERAMV1 and GDPAMV1 was found to be 8347 bp and 8011 bp, respectively (Figure 9). The mRNA folding derived from cDNA of the designed vaccines ERAMV1 and GDPAMV1 were obtained through mfold of UNAFold Web server and output provided as SS-count plot, the structure of nucleic acid, and structure dot plot ( Supplementary Figures 1  and 2). The predicted minimum free energy (DG in kcal/mol) for the top structure of mRNA of ERAMV1 and GDPAMV1 was found À685.4 and À484.90, correspondingly. At 1.0 M [Naþ] ionic condition, the obtained melting temperature for folded mRNA of ERAMV1 and GDPAMV1 was found at 78.70 C and 74.40 C, correspondingly.

Immune simulation of designed vaccines
The innate immune system plays an important role in the inhibition of fungal growth as well as in the prevention of Mucormycosis at an early stage of infection. Besides these, the adaptive immune system also provides a significant role in the clearance of fungi through the involvement of Th and Tc cells. The corresponding roles of IFN-c (Th1) and IL-17 (Th17) cytokines have been also documented in stimulating phagocytic and neutrophilic responses to provide protection against fungi Griffiths et al., 2021;Matsuzaki & Umemura, 2018). In the case of Rhizopus oryzae, IL-17 is also responsible for the release of defensin from neutrophils which regulates the inflammatory and immunological processes (Ibrahim & Kontoyiannis, 2013;Vega et al., 2011). Similarly, the T-cells of mucormycosis patients produce IL-4, IL-10, IFN-c, and IL-17 cytokines which are associated with the induction Th1 responses, limit inflammatory pathology as well as endorse long-term memory immunity, hyphal damage, and control fungal infection, correspondingly (Brown, 2011;Potenza et al., 2011). Thus, the immune simulation of designed vaccines was performed for 350 days to predict the possible induction of immune cells along with the positive control SAPN (C2). The SAPN was designed using  bioinformatics approaches and experimentally verified as a vaccine candidate (FMP014) to induce cytokines (IFN-c, TNFa, IL-4, IL-10) and antibody (IgG) in mice against Plasmodium falciparum (Kaba et al., 2018). The obtained results of C-ImmSim for both the designed vaccines indicate their potential to induce immune response after the 1st dose. There was a drastic decrease in antigen count for both the designed vaccines after 1-2 days of 1st dose along with a higher level of antibody titration (IgM, IgG1, IgG2), and Tc cell population (Figure 10) compared to the positive control (C2) (Supplementary Figure 3). The active Th concentration of ERAMV1 is slightly lower than that of control (C2) and GDPAMV1. Fortunately, the responses of dendritic cells, macrophages, and natural killer cells of both ERAMV1 and GDPAMV1 were found to be comparable to the control (C2). Likewise, the concentration of cytokines and interleukins of both designed vaccines ERAMV1 and GDPAMV1 were also found comparable to the control (C2) except the transforming growth factor-beta (TGF-b) and IL-2. The concentration of TGF-b and IL-2 against GDPAMV1 was higher than that of C2 followed by ERAMV1 (Figure 10 and Supplementary Figure  3). According to Dos Santos et al. (2021), a higher level of IFN-c and IL-2 is associated with fungal clearance and higher resistance in pulmonary mucormycosis. Further, IL-2 is also known as a marker of induction of adaptive immune response (Mizui, 2019). The important role of TGF-b and TNFa has been documented by Shao et al. (2005) during chronic fungal disease in limiting the production of chemokines and prevention in damaging tissues. TGF-b and IL-10 produced by regulatory T-cells also control the inflammatory responses and provide defense against fungus (Golubovskaya & Wu 2016). In the case of R. oryzae, IFN-c is involved in the     activation of macrophages and neutrophils, cytokine induction, hyphal damage, and protection against mucormycosis (Montaño & Voigt, 2020). Both designed vaccines ERAMV1 and GDPAMV1 showed a higher level of plasma B lymphocytes count sub-divided per isotype (IgM, IgG1, and IgG2) than that of control (C2). The concentration of active B cell population of GDPAMV1 is higher than that of control (C2) while ERAMV1 showed a slightly lower concentration than that of control (C2) (Figure 10 and Supplementary Figure 3). The above-mentioned facts correlate with the potentials for inducing humoral as well as a cellular immune responses for both the designed vaccine candidates against mucormycosis.
There was an increase in concentrations of memory B cells and active B-cells after every dose of ERAMV1 and GDPAMV1 vaccines. The elicitation of different B-cell types indicated class switching potential which reflects the longlasting immune response of ERAMV1 and GDPAMV1 vaccines. The concentrations of Th cells (active and memory) and Tc cells were also increased after each dose. The highest concentration of Th cells and Tc cells was found in the tertiary  Figure 10). Overall, when we compared to primary, secondary, and tertiary immune responses for both the designed vaccines ERAMV1 and GDPAMV1 with the experimentally validated vaccine control (C2) of Kaba et al. (2018) for 350 days, they showed significantly better results.

Discussion
Before the pandemic, mucormycosis was a threat for only immunocompromised patients such as patients of diabetes mellitus, organ transplantation, cancer, and so on. But due to the emergence of SARS-CoV-2, there is a drastic increase in mucormycosis cases and deaths (Marr et al., 2002). The world health organization (WHO) emergency report states that India is highly affected by COVID-19 associated with mucormycosis. In India, the risk of mucormycosis is 80 times greater than that of other developed countries that have an estimated 140 per million population (WHO Emergency Report, 2021). Although various antifungal drugs and surgical therapy are available for mucormycosis but the emergence of drug resistance and painful surgery procedures, there is an urgent need to develop new therapeutics and vaccines (Nagy et al., 2021). Several scientific efforts have been done to develop multi-epitope based vaccines against different fungi such as Paracoccidioides brasiliensis (Iwai et al., 2003), Candida albicans (Tarang et al., 2020), Candida glabrata (Elhasan et al., 2021), Aspergillus fumigatus (Thakur & Shankar, 2016) and so on. To date, no in silico study was reported utilizing whole proteome-based screening of R. delemar for identification of antigens and development of multi-epitope vaccines. In this regard, the immunoinformatics approach has been found cost-effective, easy, fast, specific, and reduces the risk of failure as compared to traditional methods (Iwai et al., 2003;Nanjappa & Klein 2014;Pyasi et al., 2021). Therefore, the present study explored the option to design multi-epitope-based vaccines against R.delemar using reverse vaccinology and immunoinformatics approaches. In which, the 4 GDPA were screened that contain GPI-anchoring and signal peptides which can induce an immune response as GPI-anchored proteins have been targeted for various therapeutic development including vaccine (Delgado et al., 2003;McLellan et al., 2012;Skountzou et al., 2007). Signal peptide domains have been also found with high epitope densities which can enhance the antigenicity of protein (Carmon et al., 2015;Kovjazin & Carmon, 2014).
To validate the present in silico strategy of predicting epitope ensemble, we have used 4 randomly selected ERA that can induce the immune system. The hypothetical protein RO3G_11882 is also known as CotH3 which is directly associated with host interaction and leads to host cell invasion. The CotH (CotH3) protein is highly targeted for therapeutic purposes, such as an anti-CotH antibody (Gebremariam et al., 2019). The hypothetical protein RO3G_17119 is also known as RSA protein. Further, the hypothetical protein RO3G_17119 is an antigen and has been used as a biomarker (Sato et al., 2017). The hypothetical protein G6F43_000212 (KAG1057968.1) is also known as rhizoferrin biosynthesis protein which is a fungal NRPS-independent siderophore (NIS) enzyme secreted by R. delemar (Carroll et al., 2017). The infection of mice with CTB-siderophore results in antibody elicitation and reduction in intestinal colonization of salmonella (Sassone-Corsi et al., 2016). The endopeptidase (EIE89256.1) is associated with the allergen family of fungus and can induce IgE (Sircar et al., 2012).
The TLR2 is a well-known pattern recognition receptor (PRR) expressed on the surface of host cells such as macrophages, dendritic cells, etc mainly involved in the activation of the innate immune response against fungi (Oliveira-Nascimento et al., 2012; West et al., 2006). The TLR2 also stimulated the regulatory T cells (Treg) associated with induction of IL-17 cytokine against Aspergillus fumigatus and provided protective immunity (Raijmakers et al., 2017). The hyphae of R. oryzae can induce the genes that are associated with proinflammatory functions such as NFKBIE, HMGB1, TNFA, NFB2, IL-1B, IRAK1, and IL-8 (Chamilos et al., 2008). Keeping these facts in mind, TLR2 activating mucosal adjuvant cholera toxin B subunit (CTB) was included in the designed vaccines of T-and B-cell epitopes, with suitable linkers to induce both innate as well as adaptive immune responses (Hou et al., 2014). The protein-based clinically approved strong CTB adjuvant has been used safely in human vaccines because of its potential to induce broad CD4 þ T cell responses as well as promote long-lasting protective immunity. Several approaches have been explored to couple CTB adjuvant to antigens, either by genetic fusion or by chemical manipulation, leading to strongly enhanced immune responses to the antigens (Antonio-Herrera et al., 2018;Kilander et al., 2001;Li et al., 2014;Stratmann, 2015). Thus, the present study utilized CTB adjuvant and L. lactis expression system for the design and expression of multi-epitope vaccines including linear B-cell epitopes and T-cell epitope ensemble of GDPA and ERA against mucormycosis. In the selection of the T-cell epitope ensemble following criteria were used: highest PPC value, VaxiJen score, and immunogenicity score. In which, epitopes with the same set of HLA binding alleles were removed to achieve maximum PPC value with a minimum number of epitopes. Hence, the designed vaccines ERAMV1 and GDPAMV1 contain only 23 and 16 T cell epitopes with a worldwide PPC value of 99.20% and 92.17%, respectively. This strategy of selecting epitope ensembles ultimately reduces the cost of wet-lab experimentation. Although, selection of the 16-mer length of B cell epitope in vaccine design is not universal but seems to dominate in some recent theoretical studies like multi-epitope vaccine development (Dong et al., 2020;Moming et al., 2018) and experimental studies like IgG/IgM antibody inducing B cell epitopes , Murugan et al., 2020, Pinto et al., 2021Tajiri et al., 2010) Figure 10. Result obtained by immune simulation of designed vaccines ERAMV1 and GDPAMV1. Section a, b, represents antigen profile and antibody titration; c, d, e, f, represent the B cell population; g, h represents plasma B lymphocytes count sub-divided per isotype (IgM, IgG1 and IgG2); i, j, k, l represents CD4þ T-helper lymphocytes; m, n, o, p, represents CD8þ T-cytotoxic lymphocytes count per entity-state, q, r represents natural killer cells (total count); s, t represents macrophages: total count, internalized, presenting on MHC class-II, active and resting macrophages; u, v represents Dendritic cells (DC) that can present antigenic peptides on both MHC class-I and class-II molecules; w, x represent the concentration of cytokines and interleukins titration.
PLLPGTSTTSTGPCKT (IEDB epitope ID 134472). Recently, Ras-Carmona and group have developed a BCEPS server that also included 16-mer as default peptide length for prediction of linear B cell epitope (Ras-Carmona et al., 2021). Also, the presence of discontinuous B-cell epitopes in the designed vaccines confirms the broad-spectrum antibodies responses. The innate immune responses of both the designed vaccines were evaluated through molecular docking studies with the ClusPro2.0 tool. The ClusPro is an automated web-based server for docking two interacting proteins. The server performs three computational steps as follows: (1) rigid body docking using program PIPER, (2) clustering the 1000 lowest energy docked structures using pair wise interface RMSD (IRMSD) as the distance measure, and (3) refinement of the predicted complex structures located at cluster centers by minimizing their energy. The IRMSD is defined as the RMSD of all interface backbone atoms once superposed on the native structure (interface atoms are defined as those within space is sampled on a sphere-based grid that defines a subdivision of a spherical surface in which each pixel covers the same surface area as every other pixel. The considered 70,000 rotations correspond to about 5 degrees in terms of the Euler angles. The step size of the translational grid is 1 Å, and hence the program evaluates the energy for 10 9 -10 10  (Desta et al., 2020). Analysis of docking energies revealed that the designed vaccines have more affinity toward TLR 2 compared to control C1 that potentiates the possibility of interaction with TLR2 and can induce an innate immune response. In addition, control (C1) share common interacting residues ARG-321, TYR-323, and LYS-347 of TLR2 receptor with ERAMV1 while in the case of GDPAMV1 common interacting residues were ARG-321, TYR-326, and LYS-347 (Table 6). However, the binding residue ASP-327 was found to be common in controls (CA1 and CB1) and the GDPAMV1 vaccine. The binding residue PHE-325 of TLR 2 is shared between controls (CA1) and ERAMV1. Also, the binding residue TYR-326 of TLR 2 is shared between controls (CB1 and C1) and GDPAMV1. The interaction residues of ERAMV1 and GDPAMV1 belong to non-adjuvant regions. This indicates that the epitope regions of ERAMV1 and GDPAMV1 can easily induce an innate immune system in absence of adjuvant.
Finally, the overall immune induction of designed vaccine candidates was evaluated through immune simulation webserver C-ImmSim. The online immune simulation incorporates the key facts of existing immunological knowledge, such as the diversity of specific elements, MHC restriction, clonal selection, thymic education of T cells, antigen processing and presentation, cell-cell cooperation, homeostasis of cells created by the bone marrow, hypermutation of antibodies, Figure 10. Continued maturation of the cellular and humoral response, and memory (Rapin et al., 2010). Several scientific reports have used the C-ImmSim tool for in silico immunogenicity evaluation of their multi-epitope vaccine candidates (Kar et al., 2020 andSayed et al., 2020). Therefore, the immune simulation of designed vaccines was performed for time steps of 350 days along with the control vaccine candidate (C2) of Kaba et al. (2018). The results of the immune simulation of both the designed vaccine candidates were found more effective in terms of immunological parameters counts compared to control vaccines that provide the reason for considering the induction of both innate as well as adaptive immune responses in wet-lab experimentation against ERAMV1 and GDPAMV1.
After successful evaluation of immunological parameters, the designed vaccines were evaluated for their in silico cloning and expression in the L. lactis host. Several scientific work of literature depicted the successful use of L. lactis as a host for the expression of vaccines that induce high levels of antibodies and cytokines responses against different diseases such as Helicobacter pylori (Aliramaei et al., 2019), Mycobacterium bovis (Stedman et al., 2019), Mycobacterium tuberculosis (Mancha-Agresti et al., 2017), cancer (Zhang et al., 2016), Polish avian H5H1 influenza (Szczepankowska et al., 2017) and Staphylococcus aureus (Samazan et al., 2015). Thus, the available experimental evidence supports the use of L. lactis as a safe vehicle i.e., free of endotoxins, and produces pure, properly folded, and stable proteins of fungal genes such as the cellulase gene of Neocallimastix (Ozkose et al., 2009) andribotoxin (Alvarez-Garc ıa et al., 2008). The designed vaccine candidates were successfully expressed in L. lactis as per the selected parameters of SnapGene software.

Conclusion
Mucormycosis is a neglected disease that limits the development of drugs, vaccines, and therapeutics. But due to the rapid increase in the number of cases of mucormycosis in the ongoing COVID-19 pandemic, there is an urgent need to develop an effective vaccine against mucormycosis. The present study utilized the power of reverse vaccinology and immunoinformatics approaches for the proteome-wide screening of potential antigens and multi-epitope-based designing of vaccine candidates ERAMV1 and GDPAMV1 against R. delemar in a very short duration of time. Further, the molecular docking and immune simulation of designed vaccines (ERAMV1 and GDPAMV1) provide the clue for inducing both innate and adaptive (humoral and cellular) immune responses against mucormycosis. Although maximum possible in silico characterization and evaluation of designed vaccines were performed, their experimental validation (in vitro) is required before testing into animal models (in vivo). Besides these, the strategies utilized in this study can be significantly useful for the designing of vaccines and therapeutic against other fungal or parasitic diseases.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Funding
The author(s) reported there is no funding associated with the work featured in this article.