Sequences of FTR1 protein from four different fungi were retrieved in FASTA format from the NCBI (https://www.ncbi.nlm.nih.gov/) database.

MHC class-I or CD8+ cytotoxic T-lymphocytic (CTL) and MHC class-II or helper T-lymphocytic (HTL) epitopes were targeted to get an effective vaccine. The most antigenic CTL epitopes were screened using the NetCTL server (http://www.cbs.dtu.dk/services/NetCTL/) [33]. CTL alleles were grouped into 12 superfamilies (A1, A2, A3, A24, A26, B7, B8, B27, B39, B44, B58, B62). Targeted FTR1 protein was screened against each of the CTL superfamily members [kept threshold values as HLAI binding (epitope identification) >0.75, weight on proteasomal C-terminal cleavage = 0.15, and weight on TAP (transport efficiency) = 0.05]. Higher predicted scores from the NetCTL server (>0.75) correspond to higher specificity of the epitopes (21). Next, HTL epitope binding was predicted using the IEDB recommended 2.22 prediction method of the IEDB prediction server (22) with a low percentile rank (<1.00) as the selection criteria. In this case, a lower percentile rank signifies a higher affinity of the epitopes. Finally, the human allele reference sets were utilized to predict HTL binding (23, 24). The B-cell epitope structure can be classified into two types based on the spatial structure of the epitopes: continuous (linear) and discontinuous (conformational) epitopes (25). The BepiPred linear epitope prediction method 2.0 of the IEDB prediction server (http://tools.iedb.org/bcell/) was used to identify the linear B-cell lymphocyte epitopes (LBL), keeping all parameters set to default (22). The discontinuous (conformational) epitopes, on the other hand, were predicted using the ElliPro prediction tool (http://tools.iedb.org/ellipro/) of the IEDB server (22). This method uses solvent accessibility and flexibility to predict epitopes. The ElliPro prediction tool selected antigenic residues from the predicted three-dimensional (3D) structure. The lowest score and maximum distance (Angstrom) were calibrated using the default mode with 0.5 and 6, respectively (26). As the 3D structure of FTR1 was not available, its model was generated using the RaptorX web server (http://raptorx.uchicago.edu/). The server uses an excellent and efficient template-based technique to predict query proteins’ tertiary or 3D structure (27). The GalaxyRefine server (http://galaxy.seoklab.org/cgi-bin/submit.cgi?type=REFINE) was then used to refine the modeled tertiary structure. The tool uses dynamics modeling and a CASP10-tested refinement method to produce better-refined structures (28, 29). This refined 3D model of the FTR1 protein was eventually used to predict the conformational (discontinuous) B-cell epitope. In addition, the IEDB server’s Emini surface accessibility prediction tool was utilized to predict surface epitopes from the conserved region using the 1.0 default threshold value (30).

A conservation study of the predicted epitopes was performed using the selected fungal species with an accessible annotated sequence from NCBI. More comprehensive coverage and protection against all mentioned virulent rhizopus and mucor species were achieved by selecting epitopes in the conserved regions of the R. oryzae sequence. CLC drug discovery-workbench software version 3.0 (https://digitalinsights.qiagen.com/) was used to perform multiple sequence alignments of the selected proteins to identify epitopes that are conserved across all the targeted species. The conservancy test guarantees that the polyvalent vaccine has broad-spectrum activity against the targeted species or types of fungus.

Epitopes that were highly antigenic, non-allergenic, non-toxic, 100% conserved among the desired species, and had the proper cytokine-inducing ability (for HTL epitopes) were considered the most promising epitopes for the vaccine development process. The antigenicity and allergenicity of the selected epitopes were predicted using the VaxiJen v2.0 (http://www.ddg-pharmfac.net/vaxijen/VaxiJen/VaxiJen.html) server (31) and the AllerTOP v2.0 (https://www.ddg-pharmfac.net/AllerTOP/) server, respectively (32). Afterward, the toxicity prediction of the epitopes was analyzed using the support vector machine (SVM) prediction method of the ToxinPred server (http://crdd.osdd.net/raghava/toxinpred/), with all parameters kept at default (33). The SVM method has great accuracy in distinguishing toxic and non-toxic epitopes. Finally, the TMHMM v2.0 server (http://www.cbs.dtu.dk/services/TMHMM/) was used to predict the transmembrane topology of all epitopes, with the parameters set to their default settings (34). The cytokine-inducing ability of the HTL epitopes is important as induced cytokines such as interferon-gamma (IFN-gamma), interleukin-4 (IL-4), and interleukin10 (IL-10) stimulate diverse immune cells, i.e., cytotoxic T-cells, macrophages, B-cells, and generate substantial immune responses (35). The IFN-gamma induction capability of the predicted HTL epitopes was determined using the Design Module and the Hybrid (Motif + SVM) prediction approach of the IFNepitope (http://crdd.osdd.net/raghava/ifnepitope/) server. Here, the Hybrid prediction approach is one of the highly accurate approaches for predicting the IFN-gamma inducing capability of the epitopes (36). Furthermore, the ability of the HTL epitopes to induce IL-4 and IL-10 was assessed using the IL4pred (https://webs.iiitd.edu.in/raghava/il4pred/index.php) and IL10pred (http://crdd.osdd.net/raghava/IL-10pred/) servers, respectively (37, 38). The SVM algorithm was utilized in both servers, with default threshold values of 0.2 and -0.3, respectively.

Linker peptides were used to facilitate the conjugation of the most promising epitopes to create a fusion protein. The AAY, GPGPG, and KK linkers conjugate the CTL, HTL, and BCL epitopes. The EAAAK linkers allow partitioning domains in bi-functional fusion proteins, while the GPGPG linkers can avoid junctional epitope formation and improve immune processing and presentation (39, 40). The AAY linker is also commonly employed in silico vaccine designing due to its known ability to form effective epitope conjugation (41). In addition, the bi-lysine (KK) linkers are involved in independent immunological activities of the epitopes used in vaccines (42). Additionally, a Toll-like receptor 4 (TLR4) agonist RS09 (APPHALS) was incorporated at the N-terminus of the final vaccine design. The adjuvant and linker peptides were chosen based on the study previously conducted by Pandey et al. (43). RS09 is a synthetic version of lipopolysaccharide (LPS), which acts as a natural TLR4 ligand (44). Therefore, the presence of RS09 allows for TCR co-stimulation, resulting in a more robust immune activation. Synthetic adjuvants, such as RS09 and Freund’s adjuvant, are a safer alternative and are considered an upgrade over traditional vaccination methods (45). Further, the human beta-defensin-3 was also used as an adjuvant sequence and joined the EAAAK linker sequences. The adjuvant and epitopes were also conjugated with the pan HLA-DR epitope (PADRE) sequence. The PADRE sequence boosts the immunological response by increasing the capacity of CTL epitopes of the vaccines (46). Lastly, a TAT sequence (11aa) was appended to the C-terminal of the modeled vaccine to allow intracellular delivery (47).

The antigenicity of the vaccine construct was predicted by the VaxiJen v2.0 server, keeping a threshold of 0.5 (31). Thus, the constructed vaccine should be highly antigenic to induce a fast and efficient immune response. Afterward, the AllerTop v2.0 server was used to analyze the allergenicity of the vaccine construct. The biophysical study of the built vaccine was then performed using the ProtParam service (https://web.expasy.org/protparam/) (48). MHCII-NP (http://tools.iedb.org/mhciinp/) on the IEDB server and NetChop3.1 (http://www.cbs.dtu.dk/services/NetChop/) server were used to analyze proteasome cleavage of the final BF Vaccine (BFV) construct (49, 50).

Following the biophysical analysis, secondary structure prediction of BFV was performed using the PRISPRED 4.0 prediction method of PRISPRED (http://bioinf.cs.ucl.ac.uk/psipred/) server, keeping all the parameters as default. The server efficiently predicts the proportion of amino acids in alpha-helix, beta-sheet, and coil structures (51–55). Next, the RaptorX webserver was used to predict the tertiary or 3D structure of the vaccine construct. The predicted tertiary structure of the BFV was then refined using the GalaxyRefine module of the GalaxyWEB server (http://galaxy.seoklab.org/). In addition, the SWISS-MODEL Workspace (https://swissmodel.expasy.org/assess) was used for the structural evaluation of the refined protein (56). For protein validation, another online tool, ProSA-web (https://prosa.services.came.sbg.ac.at/prosa.php), was utilized in conjunction with SWISS-MODEL workspace. The z-score generated by the server for each query protein represents the quality of the protein structures. A query protein with a z-score that falls within the range of the z-scores of all experimentally determined protein chains in the existing PDB database signifies the higher quality of the target protein (57).

The vaccine protein disulfide engineering was performed to understand the conformational stability of folded proteins, using Disulfide by Design 2 (http://cptweb.cpt.wayne.edu/DbD2/) server (58). During the analysis, the χ3 angle was set at -87° or +97° 10 while the Cα-Cβ-Sγ angle was kept at its default value of 114.6° ± 10. Residue pairs that possess energy less than 2.2 Kcal/mol were selected and altered to cysteine residues to form disulfide bridges (59). Because 90% of native disulfide bonds have an energy value of less than 2.2 Kcal/mol, 2.2 Kcal/mol was used as the threshold for selecting the disulfide bonds (58).

When TLR proteins encounter vaccines, they induce strong immunological responses to produce immunity towards that particular vaccine, thus mimicking the initial immune response against pathogenic infections (60). Therefore, the constructed vaccines should have a binding affinity with the TLRs. In protein-protein docking, the developed vaccine was docked against numerous TLRs i.e., TLR-1 (PDB ID: 6NIH), TLR-2 (PDB ID: 3A7C), and TLR-4 (PDB ID: 4G8A). Three separate online tools were used to perform protein-protein docking to increase prediction accuracy. The docking was first done with ClusPro v2.0 (https://cluspro.bu.edu/login.php), where the lower the energy score, the better the binding affinity. The energy score is calculated by the ClusPro server using the equation E = 0.40Erep + (−0.40Eatt) + 600Eelec + 1.00 EDARS (61, 62). The docking was then conducted by the PatchDock server (https://bioinfo3d.cs.tau.ac.il/PatchDock/php.php) and refined by the FireDock server (http://bioinfo3d.cs.tau.ac.il/FireDock/php.php) (63–65). Finally, the third round of docking was executed by the HawkDock server (http://cadd.zju.edu.cn/hawkdock/) in conjunction with the Molecular Mechanics/Generalized Born Surface Area (MM-GBSA) study (66). Lower output scores in these servers translate to higher binding affinity and vice versa (67). Discovery Studio Visualizer was used to visualize the most effective vaccine-TLR combination (68).

The conformational B-cell epitopes of the constructed vaccine were predicted using the IEDB ElliPro tool (http://tools.iedb.org/ellipro/), with all parameters set to default (22). In addition to cell-mediated immunity, humoral immunity is essential in fighting against infections inside the human body. The body’s humoral immunity depends on B-cells, which secrete antibodies in response to antigens. As a result, the vaccine construct should have effective conformational B-cell epitopes to boost humoral immunity.

Molecular Dynamic (MD) simulations were performed using the GROMACS 2020.3 simulation package (69). The AMBER99SB-ILDN force field was combined with the TIP3P water model (70). Each system received 0.15 M KCl. Potential energy minimization with a step size of 1 fs to the maximum force of 1000.0 kJ/mol/nm was used to relax the structure and avoid steric clashes in subsequent simulations. The system’s pressure and temperature were equilibrated to 1 atm and 310K, respectively, by running NVT and NPT simulations (100-ps each). The system’s pressure and temperature were regulated using a modified Berendsen thermostat and a Parinello-Rahman barostat with time constants tau t = 0.1 ps and tau p = 2 ps, respectively (71, 72). Productive 200 ns MD simulations with a 2-fs time-step were performed in the isothermal-isobaric ensemble for all six systems. The LINCS algorithm constrains the hydrogen-atom bonds (73). The Particle-Mesh Ewald summation scheme was used to calculate long-range electrostatic interactions (74).

Using the C-ImmSim server (https://kraken.iac.rm.cnr.it/C-IMMSIM/index.php), an immune simulation analysis of the vaccine construct was conducted to estimate its immunogenicity and immune response profile. The server uses machine learning algorithms and a position-specific scoring matrix (PSSM) to predict real-life-like immunological interactions (75). All parameters except the time steps were left at their default levels during the analysis. The number of simulation steps was set to 1050, while the time steps were kept at 1, 84, and 170. Because the recommended period between two subsequent doses of most commercial vaccinations has been proven to be four weeks, administration of three injections would need to be spaced four weeks apart (76). The figures were used to calculate the Simpson’s Diversity Index, D.

Because an amino acid might be coded by more than one codon in different organisms, codon adaptation is used to predict the most appropriate codon that efficiently encodes the required amino acid in a particular organism. The Java Codon Adaptation tool (http://www.jcat.de/) was used to optimize the codon adaptation study of the vaccine protein, and the prokaryotic Escherichia coli (E. coli) strain K12 was chosen as the target organism for increasing the expression efficiency of the final vaccine protein (77). In the server, rho-independent transcription termination, prokaryotic ribosome binding site, and cleavage site of restriction enzymes were avoided to ensure the vaccine gene’s accurate translation. Instead, XhoI and BamHI restriction endonuclease sites were added to the vaccine’s N and C terminals, respectively. The modified DNA of the constructed vaccine was then incorporated between the XhoI and BamHI genes in the pET28a (+) vector using the SnapGene restriction cloning software (78). The ubiquitin-like modifier or SUMO tag and 6XHis-tag in the pETite vector plasmid aid in the recombinant protein’s solubilization and viable affinity purification (79).

Two online tools, Mfold version 2.3 (http://www.unafold.org/mfold/applications/rna-folding-form-v2.php) and RNAfold (http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi) were used to predict the mRNA secondary structure prediction of the vaccine gene. Both of these servers thermodynamically predict mRNA secondary structures and report the least free energy (ΔG Kcal/mol) for each of the structures created. The folded mRNA is more stable when the minimum free energy is lower and vice versa (80–83). To analyze mRNA folding and secondary structure of the vaccine, the DNA<->RNA->Protein tool (http://biomodel.uah.es/en/lab/cybertory/analysis/trans.htm) was used to convert the newly adapted DNA sequences from the JCAt server into its probable RNA sequence. The generated RNA sequence was then copied and pasted into the Mfold and RNAfold servers for prediction, with all parameters set to default.

FTR1 protein of four different virulent funguses was selected after reviewing the literature from the NCBI database (Supplementary Table S1). The selected fungi were Rhizopus oryzae, Rhizopus delemar RA 99-880 (a reclassified Rhizopus oryzae), Rhizopus azygosporus, Rhizopus stolonifer, and Mucor circinelloides 1006PhL. The FTR1 protein sequence of the Rhizopus delemar RA 99-880 was the only reviewed protein sequence in the Uniprot database (Entry name: FTR1_RHIO9) (https://www.uniprot.org/), considered as the model sequence to carry out the experiments.

For epitope mapping, the sequence of the FTR1 protein of R. delemar was evaluated, and a total of 309 epitopes were screened out for further research. After screening the FTR1 protein against all the CTL superfamilies in NetCTL-1.2 Server, epitopes were predicted based on the predicted scores. A total of 182 epitopes with higher scores (>1.00 signifies stronger binding affinity) were enlisted (Supplementary Table S2). Further, the top 122 HTL epitopes with low percentile scores (<1.00 represents better binding affinity) in the IEDB server were chosen (Supplementary Table S3). Finally, the LBL epitopes were predicted using the BepiPred linear epitope prediction method 2.0 of the IEDB prediction server. The binder’s protein score to B-cell was an average of 0.450, a minimum of 0.189, and a maximum of 0.693. Since all values were equal or higher than the default threshold of 0.500, all the epitopes were considered LBL epitopes (Supplementary Figure S1). Additionally, the average surface accessibility area of the FTR1 protein was found to be 1.000, a maximum of 5.900, and a minimum of 0.025 in Emini’s surface accessibility prediction test for a potent LBL epitope. Values ≥ the default threshold of 1.000 were selected (Supplementary Table S3). In the end, 5 LBL epitopes of 8 to 20 AA were chosen for further analysis (Supplementary Table S4). The ElliPro prediction server filters out antigenic residues from the 3D model of the FTR1 protein to predict the conformation of B-cell epitopes (Table 1). The lowest score and maximum distance (Angstrom) in the default mode were calibrated to 0.5 and 6, respectively.

Fully conserved epitopes of the targeted fungal species were further used for the conservation analysis. Afterward, antigenicity, allergenicity, toxicity, cytokine inducing ability, and transmembrane topology of the identified epitopes were analyzed. Epitopes found to be highly antigenic (>0.5), non-allergenic, non-toxic, 100% conserved among selected species, and strong cytokine inducing ability (for HTL epitopes) was considered the most promising epitopes for the vaccine development process (Table 2).

In addition, two HTL epitopes were found to have low percentile scores (less than 1.00) and appreciable cytokine-producing capability (IFN-gamma, IL-4, and IL-10) (Table 3).

The Emini score was considered for selecting LBL epitopes, which revealed KTERMQEKWKVK as the best sequence with a score of 3.077, highly antigenic, non-allergenic, and non-toxic. Therefore, we selected eleven CTL, HTL, and LBL epitopes as the most promising epitopes for the peptide fusion, fully conserved among the selected species (Figure 3).

The vaccine was constructed using the eleven most promising epitopes. RS09 and beta-defensin-3, appended at the N-terminal of the vaccine construct, were used as an adjuvant. In the C-terminal of the vaccine, the PADRE and TAT sequence were conjugated to act as strong immunity inducers. EAAAK, AAY, GPGPG, and KK linkers were employed in the correct positions to connect the epitopes. The newly developed vaccine is designated as the BF vaccine (BFV). Figure 4 depicts a schematic representation of the constructed vaccine BFV.

The newly constructed BFV was potent antigen and non-allergenic (Supplementary Table S5). The biophysical studies revealed that BFV possesses a high (primary) theoretical pI (9.97). The chemical formula of the vaccine protein can be denoted as C₁₃₂₃H₂₀₃₄N₃₄₂O₃₂₂S₁₀, contains 4031 atoms, and has a molecular weight of 28203.40. It has a half-life of 4.4 hours in mammalian reticulocytes in vitro, >20 hours in yeast in vivo, and more than 10 hours in the E. coli cell culture system, which was satisfactory. Moreover, the GRAVY value of the vaccine design was 0.111, indicating that the protein is hydrophobic. The extinction coefficient was 80705 M^-1 cm^-1 at 280 nm measured in water. BFV was also shown to have a high solubility score of 0.640 in the Protein-sol server. Proteasome cleavage analysis using MHCII-NP from the IEDB server and NetChop3.1 server revealed that BFV was found to be cleaved by proteasomes to give rise to the predicted T- cell epitopes (Supplementary Tables S6, S7).

The secondary structure prediction of BFV displayed by the PSIPRED server revealed that its alpha-helix structure comprised the highest percentage of amino acids while the lowest proportion was discovered in the β-strand (Supplementary Figure S3). The RaptorX server was used to generate the 3D structure of BFV, which was then refined so that the 3D BFV protein could closely resemble the native protein structures. The refined 3D structure of BFV was further evaluated by the Ramachandran plot generated by the SWISS-MODEL Workspace and validated by the ProSA-web server. The Ramachandran plot analysis showed that BFV had a 96.43% favored region with 0.00% outliers and included A254 LYS A23 ARG in the Rotamer outliers (1.09%). In addition, BFV has a z-score of -4.38, which is within the range of all experimentally established X-ray crystal structures of proteins in the Protein Data Bank. Moreover, protein validation investigation confirmed the high quality of the revised structure (Figure 5).

The Disulfide predicted the disulfide bonds of BFV by Design 2 server. Based on a selection criterion, the server determines the possible sites within a protein structure with the greater possibility of pairs of amino acids that can form disulfide bonds. Only amino acid pairs with a bond energy of less than 2.2 kcal/mol were chosen for this experiment. Four amino acid pairs in BFV protein were found to have bond energy less than 2.2 kcal/mol: 23 Arg and 59 Ala, 62 Ala and 65 Val, 119 Trp and 122 Val, and 213 Lys 216 Arg. These selected amino acid pairs were used to form the mutant version of the original BFV by forming disulfide bonds in the Disulfide by Design 2 server (Supplementary Figure S4).

The protein-protein docking study was performed to evaluate the ability of BFV to interact with different TLRs. Three online tools were used to perform protein-protein docking to improve prediction accuracy. When docked by ClusPro 2.0, it demonstrated very high binding affinities with all of its targets. In PatchDock and HawkDock, BFV has also shown excellent interaction with all of its targets. As a result, it may be argued that the BFV should elicit a strong and consistent immune response when administered (Supplementary Table S8). Finally, the discovery Studio Visualizer shows the interactions of the BFV with several TLRs (Figure 6).

The conformational B-cell epitope prediction of BFV revealed seven regions within the vaccine with scores ranging from 0.566 to 0.913, spanning a total of 120 amino acid residues (Figure 7). Supplementary Table S9 shows the individual scores of each of the conformational B-cell epitopes.

For 200 ns simulations, trajectory analysis was performed. Since the periodic boundary conditions were met during the simulation, the protein molecules in the complex were re-centered and returned to the simulation cell. Following that, an initial analysis of the trajectories was performed, which included calculating the RMSD, radius of gyration, and RMSF. Figure 8 depicts the corresponding graphs. Figure 8A depicts the RMSD changes. Figure 8A shows that the RMSD stabilizes around 0.8 for BF1, around 1.0 for BF2, and around 0.6 for BF4. Given the size of the simulated systems, this indicates that they are very stable. Figure 8B shows that the radius of gyration for each complex tends to be stable, indicating that the compactness of complexes does not change significantly during simulation. Figure 8C depicts the RMS fluctuations for the complexes’ C-alpha atoms. According to the graphs, the mobility of the receptor atoms is significantly lower than that of the vaccine atoms.

The BFV elicited both robust primary and secondary immune responses. After administration of the three doses of BFV, increased levels of active B cells (Figures 9B, C), plasma B-cell (Figure 9D), helper T-cell (Figure 9E, F), and cytotoxic T-cell (Figures 9G, H) along with gradually elevated levels of different immunoglobulins (Figure 9A) were predicted. At the same time, the vaccine stimulates helper T-cells which eventually result in an improved adaptive immunity (84, 85). Again, improved immunological memory development, greater antigen clearance rate, increased dendritic cells, and macrophages suggest outstanding antigen presentation by the Antigen Presenting Cells, i.e., dendritic cells and macrophages (Figures 9I, J). Moreover, the vaccines induce diverse cytokines such as IFN-gamma, IL-23 (interleukin-23), IL-10, and IFN-beta (interferon-beta), that instigate immune response against pathogens and defend the body (86, 87). BFV was also shown to induce different types of cytokines (Figure 9K). Simpson’s Index (D) was measured as a nominal value, signifying that the vaccine has a more diverse effect in the analysis (88). Thus, the immune simulation study findings strongly suggest that BFV is a promising vaccine that stimulates a robust immunogenic response against FTR1 protein after its administration.

The protein sequence of BFV was adapted in the JCat server for codon adaptation. BFV had a codon adaptation index (CAI) value of 0.798, suggesting that the DNA sequences of BFV contain a higher proportion of the codons that are expected to be used in the cellular machinery of the target organism E. coli strain K12 (Supplementary Figure S5). As a result, it was predicted that the E. coli strain K12 would produce the BFV efficiently (80, 81). Furthermore, the produced sequence had a GC content of 49.74%. Afterward, the predicted DNA sequence of BFV was introduced into the pETite vector plasmid between the XhoI and BamHI restriction sites. Figure 10 depicts the newly created recombinant plasmid. SUMO tag and 6X His tag of the vector plasmid facilitate vaccine purification during the downstream processing (89). Following that, the Mfold and RNAfold servers predicted the secondary structure of the BFV mRNA. The Mfold server calculated a minimum free energy score of 199.10 kcal/mol, consistent with the RNAfold server’s -225.80 kcal/mol prediction. Supplementary Figure S6 depicts the secondary structure of BFV mRNA.