Introduction

Bacillus licheniformis is a Gram-positive spore-forming bacterium, with a wide range of beneficial and pathogenic properties characterizing its history. It is acknowledged as a probiotic strain that can strengthen human health through its consumption in limited amounts.¹ However, there have been documented instances of B. licheniformis inducing human and bovine infections.^2,3 It is a common contaminant of food products and can cause spoilage and foodborne illness.⁴ In addition, B. licheniformis is implicated in infections in immunocompromised patients, potentially threatening public health.⁵ The significance of One Health approach is highlighted by the emergence and pathogenicity of B. licheniformis, including the interrelationships in humans, animals, and the environment. Infections caused by B. licheniformis not only entail a potential hazard to human health but also impact livestock populations and ecological systems.⁶ The dissemination of this bacterium can transpire via human–animal interaction, in addition to the pollution of environmental reservoirs, such as water and soil.

Despite the availability of conventional antimicrobial therapies and probiotic combinations^2-3,6 to manage B. licheniformis infections, these strategies face significant limitations. Antimicrobial therapies, while effective in some cases,⁷ contribute to the emergence of antimicrobial resistance, reducing their long-term utility. Similarly, probiotic combinations aimed at minimizing pathogenicity⁸ are not universally effective and may pose risks of infection in immunocompromised patients. Current interventions⁹ also fail to address the bacterium’s ability to persist in environmental reservoirs, complicating its eradication.⁶ These limitations highlight the need for novel targeted approaches to mitigate the risks posed by B. licheniformis.

The economic implications of the pathogenicity of B. licheniformis are also important to discuss. According to estimates, direct expenses associated to manage B. licheniformis infections can be significant, with an average cost of $10,000 per treatment in the United States.¹⁰ Moreover, indirect expenses also exist linked with these infections, such as reduced productivity and income, as those infected may be incapable of attending school or work. In addition, the image and standing of probiotic merchandise and corporations may suffer adverse consequences because of B. licheniformis infections, resulting in declined revenue and earnings.⁸ These economic ramifications highlight the significance of tackling the pathogenicity of B. licheniformis through various measures, including creating an in-silico vaccine. Such interventions can alleviate financial strains and foster human health and economic prosperity.

Epitope-based vaccines offer unique advantages over conventional vaccine development strategies. By selectively targeting immunogenic and nonallergenic protein regions, they allow for a high degree of specificity, reducing the risk of adverse immune reactions.^11,12 This precision is particularly important for a bacterium such as B. licheniformis, whose dual role of being a probiotic strain necessitates careful balancing of its beneficial and pathogenic characteristics. Additionally, in-silico approaches accelerate the identification and evaluation of vaccine candidates, enabling the designing of safer and more effective interventions.¹³ These attributes make epitope-based vaccines a promising solution to address the multifaceted challenges posed by B. licheniformis infections.

The present investigation aims to surpass the limitations encountered in the conventional vaccine development process for B. licheniformis by using a computer-simulated approach. We have selectively targeted nonallergenic but antigenic proteins of B. licheniformis through immunoinformatic analysis. This approach presents numerous benefits, such as the capacity to evaluate a wide range of prospective vaccine candidates and anticipates their immunogenicity without requiring extensive laboratory experimentation.¹⁴ In addition, in silico design enhances vaccine effectiveness, safety, and stability, resulting in a more precise and individualized strategy for addressing B. licheniformis infections.¹⁵

This study addresses the challenges posed by B. licheniformis as probiotic and foodborne pathogen by presenting a targeted vaccine strategy that ensures safety without compromising efficacy. By integrating 16S rRNA-based molecular identification with immunoinformatic screening of antigenic nonallergenic epitopes, the work establishes a rational framework for vaccine design that overcomes the limitations of conventional approaches. This contribution advances computational vaccinology by demonstrating how in silico methods can accelerate candidate selection, enhance safety profiles, and reduce dependency on extensive laboratory testing. Importantly, the findings extend beyond methodological innovation, aligning with the One Health perspective by offering a preventive solution that addresses human, animal, and environmental dimensions of B. licheniformis infections.

Materials and Methods

A fusion of in vitro bioinformatic and immunoinformatic approaches was used to construct a vaccine candidate. These approaches and thresholds were chosen based on the previously published studies.^16-18 The methodology is outlined in Figure 1.

Figure 1 Roadmap to constructing a vaccine candidate against B. licheniformis.

Isolation and identification of bacterial strain

The bacterial strain under investigation was isolated from soil samples collected from the rhizosphere of Psidium guajava (guava), Polyalthia longifolia (false or ulta Ashoka), and kimchi with Brassica oleraceavar, capitata (cabbage) being the main ingredient. The cultured strains were led to biochemical tests, including Gram staining, catalase, gelatin, sugar fermentation, and 6.5% NaCl to identify their lineage. Genomic DNA was extracted from bacterial cells using the detergent cetyltrimethylammonium bromide (CTAB) method, and the 16S rRNA gene was amplified using universal prokaryotic primers 27F and U1492R.^19,20 The polymerase chain reaction (PCR) reaction consisted of 35 cycles with specific temperature profiles shown in Supplementary Table S1. The amplified products were visualized on 1% agarose gel, and sequencing was performed by 1st BASE: (https://base-asia.com/) (Selangor, Malaysia). The strain (with GENBANK accession number: OL477474.1) was identified and validated through NCBI BLAST, Ez Taxon, and MEGA-X phylogenetic analysis.

Pan-genome analysis

The study of pan-genomes involves comprehensive analysis of the complete set of genes present in each group of organisms. This approach allows for a more thorough understanding of genetic diversity and evolutionary relationships within and between species. Genomic sequences of 50 B. licheniformis strains were acquired from the GenBank repositories of the National Center for Biotechnology Information (NCBI) database. The bacterial pan-genome analysis (BPGA) was conducted through USEARCH (a bioinformatics tool) configuration for pan-genome analysis. The study involved pre-processing stages, wherein 39 strains were filtered for subsequent analyses. Sequence data were then generated to establish a sequence identity, using a threshold score of 50% (selected to balance sensitivity and specificity). The compiled results were utilized to ascertain the frequency of equivocal genes and novel gene clusters, and to compute the overview of pan-genome. The BPGA tool was utilized to obtain main sequence, which was then processed for clustering through the CD-HIT web server (http://weizhong-lab.ucsd.edu/cdhit_suite/) using a percentage identity cut off 0.5%, following the protocol established by Fu et al.²¹

Target proteins identification

The bacterial strain’s entire proteome was obtained from the NCBI genome database (https://www.ncbi.nlm.nih.gov/datasets/taxonomy/1402/). The core sequence derived from the pan-genome assessment underwent filtration and clustering procedures to eliminate redundancy. This was accomplished through the utilization of CD-HIT web server (http://weizhong-lab.ucsd.edu/cdhit_suite/). The nonredundant protein dataset was searched against the core virulence factor database (VFDB) using BLASTp (https://blast.ncbi.nlm.nih.gov/Blast.cgi?PAGE=Proteins), selecting proteins with sequence identity <30% and bit scores >100 (thresholds ensure specificity to virulence factors). The localization of the proteins under investigation within the cell was determined by means of PSORTb 3.0 (https://www.psort.org/psortb/) and subsequently verified using the CELLO2GO (http://cello.life.nctu.edu.tw/cello2go/) localization predictor tools, following the guidelines provided by Yu et al.^22,23, respectively. We further reviewed the literature²⁴ for B. licheniformis toxin proteins, cross checked their localization and verified it in the finalized core sequence. Potential biases were minimized by using multiple localization predictors and literature validation. The final protein candidates were selected because they are virulence-associated, surface-exposed, or extracellular (thus accessible to host immunity), and non-homologous to human proteome and microbiome, minimizing cross-reactivity. These combined features ensured that only relevant, safe, and immunogenic proteins were advanced for epitope prediction and vaccine construct design.

Epitope prediction

The Immune Epitope Database and Analysis Resource (IEDB) server was utilized to predict the epitopes of B-cells and T-cells for the purpose of vaccine construct design, following the modules of Vita et al.²⁵ The linear B-cells epitope prediction search tool can be accessed (http://tools.iedb.org/main/bcell/). The method of epitope prediction was employed to assess potential B-cell epitopes. The T-cell epitope prediction was carried out using the major histocompatibility complex (MHC)-I epitope prediction tool, which was accessed through at http://tools.iedb.org/mhci/. The MHC-II epitope prediction tool (accessible at http://tools.iedb.org/mhcii/) was utilized to forecast MHC-II epitopes derived from the chosen proteins. NetMHCpan server was employed for MHC class I epitope prediction, while NetMHCIIpan was used for MHC class II epitope prediction. These tools utilize advanced artificial neural networks and position-specific scoring matrices to predict peptide-MHC binding affinity with high accuracy. Epitopes of MHC-I and MHC-II were chosen for further analysis based on their half maximal inhibitory concentration (IC50) values, which were less than 50 nM (threshold ensures high binding affinity). Potential biases, such as inaccuracies in prediction algorithms, were minimized through reliance on experimentally validated thresholds and cross-comparison of predicted epitopes across tools.

Immunocompatibility of epitopes

The process of epitope selection is a critical aspect of vaccine design, wherein the epitopes must possess distinct attributes, such as nonallergenicity, nontoxicity, non-human homology, and antigenicity. To predict the antigenic potential of proteins, the Vaxijen Tool²⁶ (accessible at http://www.ddg-pharmfac.net/vaxijen) was employed by keeping the threshold value of 0.5. Additionally, the allergenic capacity was assessed using the bioinformatics-based tool called AllerTop²⁷ (available at https://www.ddg-pharmfac.net/AllerTOP/). The ToxinPred tool²⁸ (at http://crdd.osdd.net/raghava/toxinpred/) evaluated toxic epitopes. Only nonallergenic, nontoxic, non-human homologous, and antigenic epitopes were selected. This approach effectively mitigates the potential risks of adverse reactions and allergic responses. In addition, epitopes not homologous to humans can increase specificity by guiding immune responses toward the intended pathogen or bacterial strain.²⁹ Lastly, antigenic epitopes can elicit a prompt immune response, leading to the generation of antibodies and activation of immune cells. By considering these attributes, researchers can formulate safe, precise, and capable of inducing immune responses, thereby enhancing their effectiveness in eliciting protective immunity. Thresholds for each tool were set based on established literature to ensure reliability and safety.

Population coverage analysis

The IEDB population coverage tool (available at http://tools.iedb.org/population/) was employed to evaluate the proportion of individuals who are expected to produce a reaction against the chosen epitope sets. The tool considers factors such as human leukocyte antigen (HLA) genotype, T-cell restriction data, and MHC interaction, as described by Bui et al.³⁰ A global analysis of population coverage was performed. Thresholds for evaluating population coverage were based on cumulative coverage values provided by the tool, ensuring that selected epitopes target a significant portion of the global population. Potential biases, such as overrepresentation or underrepresentation of specific HLA genotypes in the database, were mitigated by including global datasets and avoiding geographic or ethnic restrictions in the analysis. This approach ensures a broader applicability of the vaccine construct, minimizing the risk of excluding any demographic group.

Design of vaccine candidate

The epitopes were integrated using linker sequences designed to elicit the intended immune response. The guidelines of Chen et al.³¹ were followed to select and design linkers, ensuring structural integrity and immunogenicity of the vaccine. Linker sequences were chosen based on the inclusion of a disulfide cleavage site, a threshold that promotes the efficient release and presentation of epitopes. The final vaccine construct was supplemented with a Pan-DR epitope (PADRE) sequences and helper sequences to enhance stability and immunogenicity. Potential biases in linker selection, such as sequence-specific immune responses or unintended interactions, were minimized by adhering to well-established design principles and experimental validation from literature.^17,32

Immunogenicity of vaccine candidate

After acquiring the vaccine sequence derived from B. licheniformis epitopes, its immunogenicity was evaluated prior to further development. The VaxiJen server²⁶ was used to predict antigenicity, with a threshold score of 0.4 set to ensure a strong antigenic response while minimizing false positives. The AllerTOP v2.0 tool²⁷ at https://www.ddg-pharmfac.net/AllerTOP/ assessed the likelihood of allergenicity, and the TOXINPRED server²⁸ at http://crdd.osdd.net/raghava/toxinpred/ evaluated toxicity.

Assessment of physicochemical properties

The physicochemical properties of the vaccine candidate sequence derived from B. licheniformis epitopes were analyzed using PROTPARAM (https://web.expasy.org/protparam/).³³ Key parameters, including half-life, stability index, molecular weight, theoretical isoelectric point (pI), and Grand Average of Hydrophathicity (GRAVY), were computed to assess the immunogenicity and stability of vaccine candidate. Thresholds for stability index (<40 for stable proteins) and GRAVY (negative values indicating hydrophilicity) were considered to ensure appropriate vaccine characteristics. The SOLpro tool from the SCRATCH protein predictor (https://scratch.proteomics.ics.uci.edu/) was used to predict disordered regions and solubility, as described by Magnan et al.³⁴. Potential biases were mitigated by combining multiple physicochemical parameters and aligning predictions with the established experimental data.

Homology and human uptake prediction

In order to evaluate the druggability and homology of the potential vaccine candidate derived from B. licheniformis epitopes, the PBIT server (http://www.pbit.bicnirrh.res.in/) was utilized as described by Shende et al.³⁵ The server assessed the potential tolerance of the host immune system and gut microbiota toward vaccine candidate. Homology analyses were conducted against human, mice, and gut microbiota sequences, with a threshold for homology set below 30% to minimize off-target immune responses and ensure specificity. Potential biases were mitigated by validating homology thresholds against established datasets and including multiple host models to ensure comprehensive analysis.

Secondary structure prediction

We utilized PSIPRED protein analysis server (available at http://bioinf.cs.ucl.ac.uk/psipred/) to predict the secondary sequence structure of vaccine candidate. The study conducted by McGuffin et al.³⁶ yielded significant findings regarding the comprehensive organization of the server and its applications. These findings included insights into disordered regions, transmembrane (TM) helix packing, fold recognition, contact analysis, and modeling of the structure’s functions and domains.³⁶ The utilization of this tool facilitated the clarification of the precise configuration of strands, helices, and coils present in the vaccine protein, thereby augmenting the all-encompassing comprehension of its structural attributes.

Immune simulation

The vaccine candidate’s potential to elicit immune responses in the human immune system was predicted using the C-IMMSIM server (accessible at https://kraken.iac.rm.cnr.it/C-IMMSIM/index).³⁷ Default parameters were maintained throughout the analysis, including the injection of vaccine candidate without lipopolysaccharide (LPS), to simulate a non-endotoxin condition. Thresholds of immune parameters, such as cytokine levels and immune cell durability, were evaluated to determine the vaccine’s potential efficacy. Potential biases were mitigated by maintaining consistency in input parameters and validating the results against known immune response patterns.

Tertiary structure prediction and validation

The 3Dpro server, SCRATCH (https://scratch.proteomics.ics.uci.edu/explanation.html#3Dpro) was utilized to determine vaccine candidate’s tertiary structure, as reported by Cheng et al.³⁸ The model exhibiting the highest confidence score was chosen from the generated models for subsequent refinement. The researchers used the GalaxyRefine tool, accessible on GalaxyWEB server (https://galaxy.seoklab.org/), to enhance the model. The tool was employed by considering various parameters, including optimal Ramachandran (RAMA) score, minimal Root Mean Square Deviation (RMSD) score, and suitable MolProbity score, as described by Ko et al.³⁹ The selection of a refined model was conducted with great care to ensure its structural integrity and quality, thereby yielding a dependable depiction of the tertiary structure of vaccine candidate.

Epitope validation analysis

The ELLIPRO tool (http://tools.iedb.org/ELLIPRO/) was used to predict antibody epitopes, providing insights into regions likely to elicit an immune response. The tool identifies conformational and linear epitopes based on structure and antigenicity scores, with a threshold score of 0.5 used to ensure specificity. Additionally, the NETCHOP tool (http://tools.iedb.org/netchop/) was utilized to predict proteasomal cleavage sites, aiding in the identification of potential antigenic fragments, as described by Kim et al.⁴⁰

Identifying molecular targets for vaccine candidate

This research identified molecular targets for vaccine candidates using receptors integral to the immune response. The Protein Data Bank (PDB) (available at https://www.rcsb.org) provided the molecular targets used for analysis in this investigation. The chosen targets included the Toll-like receptor, TLR1 (PDB ID: 6NIH), TLR2 (PDB ID: 1FYW), TLR5 (PDB ID: 3J0A), CD8 receptor (PDB ID: 1CD8), CD4 (PDB ID: IWIP), MHC-I and MHC-II (PDB ID: 2XPG and 6T3Y, respectively), MHC-I, and MHC-II. Through this analysis, we intended to learn about the vaccine candidate’s possible interactions and structural features with these proteins.

Molecular docking

Hdock (http://hdock.phys.hust.edu.cn/), a computational tool for protein–protein docking, was used to analyze the interaction between finalized epitopes and their corresponding MHC receptors, as described by Remmert et al.⁴¹. Epitopes interacting with B-cells, MHC-I, and MHC-II were docked with their respective MHC receptors to evaluate binding affinity and interaction sites. A docking score threshold of ≤-7.0 kcal/mol was applied to identify high-affinity interactions. Subsequently, epitope–MHC complexes were docked with CD8 and CD4 receptors to investigate the binding patterns of MHC molecules with T-cells. The vaccine construct was also docked with three human TLRs to identify bacterial pathogen-associated molecular patterns (PAMPs). Potential biases were addressed by cross validating docking results with known structural data and refining input parameters to enhance prediction reliability.

Molecular dynamics (MD) simulations

Molecular dynamics simulations were performed using the iMods online server (https://imods.iqfr.csic.es/) to evaluate the flexibility and stability of docking complexes, as described by López-Blanco et al.⁴² TLR-vaccine complexes and other relevant structures were subjected to computational analysis to predict parameters, such as the energy required for complex denaturation and RMSD values, with a stability threshold of RMSD ≤ 2.0 Å considered optimal. This analysis provided critical insights into conformational stability and interaction dynamics of the complexes, aiding in the optimization of vaccine candidate’s design.

Results

Isolation and identification of B. licheniformis

The isolated samples and Petri plates, depicted in Figures 2A and 2B, demonstrated the efficacy of isolation and purification process. All biochemical tests were positive, indicating that the isolated strain belonged to the Bacillus genus. The successful 16S rRNA amplification was confirmed through visualization of amplified products on a 1% agarose gel, as illustrated in Figure 2C. The sequencing of 16S rRNA plays a crucial role in identifying the bacterial species present in each sample, which helps in pinpointing the pathogens that may be targeted by a vaccine. By identifying the specific strains of bacteria involved in infections, we can focus on those that pose maximum threat to public health. This enables to take informed decisions about which strains to include in vaccine formulations. Additionally, post-vaccination, 16S rRNA sequencing can help monitor bacterial populations, ensuring the vaccine remains effective against circulating strains. Thus, the identification of pathogens through 16S rRNA sequencing directly informs the vaccine development process, from pathogen selection to monitoring vaccine efficacy. Following this, the amplified products underwent sequencing analysis conducted by 1st BASE (Selangor, Malaysia). Subsequent examination utilizing the NCBI BLAST repository, and the Ez-Taxon computational resource evidenced that the sequence manifested a noteworthy resemblance with various B. licheniformis strains. The BLAST findings revealed a considerable resemblance between both acquired sequence and established B. licheniformis sequence, implying that most bacterial isolates were affiliated with this species in the Bacillus genus. The phylogenetic tree, shown in Figure 2D, depicted evolutionary interrelationships among 15 bacterial isolates. A bootstrap value of 1000 was employed to assess confidence of the phylogenetic tree, which provided substantial evidence for positioning of isolates within the Bacillus genus. The phylogenetic tree visually represented genetic relatedness and clustering patterns among isolates.

Figure 2 Isolation and identification of the strain. (A) Isolated samples; (B) Petri plates containing selected purified cultures; (C) 16S rRNA amplification gel image; and (D) phylogenetic tree.

Pan-genome analysis

The core proteomes of 39 strains of B. licheniformis, undergone sequencing, collectively identified 93,678 proteins. Figure 3 presents the pan-genome analyses of the bacterium. The core proteomic dataset was subjected to a redundancy filter offered by the CDHit server, resulting in the identification of 1476 distinct proteins, each represented individually. The nonredundant proteins were subjected to further analysis for subcellular localization, which resulted in the identification of 32 lipoproteins, 195 membrane proteins, 75 secreted proteins, and 1168 cytoplasmic proteins.

Figure 3 BPGA analysis of B. licheniformis. The figure represents the (A) number of new genes; (B) distribution of gene families; (C) core-pan plot; and (D) KEGG distribution.

Target proteins identification

Of the 75 proteins secreted, 49 were identified as having the potential to elicit an immune response, satisfying the requirements for being regarded as potential targets for vaccine development. Subsequent examination revealed that 38 proteins with pathogenic properties exhibited a molecular weight <110 kDa, with 22 possessing either one or no TM helices. This observation highlights the necessity for further inquiry.

Immunogenicity of the 22 proteins that possess one or fewer TM helices was subjected to further assessment. The VaxiJen server was employed to determine the antigenicity of these proteins, and it identified 12 of them as being antigenic. It was determined that eight proteins exhibited nonallergenic and adhesive properties, rendering them viable candidates for vaccine targets. Significantly, it was observed that all of the proteins exhibited no homology with the human proteome. Thorough examination of literature determined that two of the eight proteins (discussed in Table S2) met the criteria to be considered as potential vaccine candidates.

Epitope prediction

The Immune Epitope Database (IEDB) generated eight B-cell epitopes from both proteins. Water solubility, antigenic, and flexible nature of identified five B-cell epitopes were characterized. Furthermore, a substantial quantity of epitopes was predicted for MHC class I and MHC class II, comprising roughly 45,000 and 50,000 epitopes. A preliminary analysis was conducted to select a subset of 20 epitopes for each MHC class with an IC₅₀ value >50 nM. Following this, a screening procedure was conducted, considering the antigenicity, lack of allergenicity, absence of toxicity, and dissimilarity to human for every epitope. Consequently, nine epitopes were selected as ultimate targets from each protein. Table S3 comprehensively summarizes the epitopes and their corresponding target MHC alleles.

Population coverage analysis

The epitopes’ collective coverage was predicted to be 89.42% of the populace. Notably, the assessment of MHC class II allele coverage was restricted due to the employed tool’s incomplete forecast of most alleles. Hence, the real population coverage of MHC class II alleles is anticipated to surpass the current prediction.

Design of vaccine candidate

The sequence of vaccine candidate (shown in Figure 4A) adhered to a specific design, incorporating diverse elements. The sequence was initiated with an EAAAK linker, which functioned as a spacer region. It also augments the stability and immunogenicity of the construct. Integrating B-cell epitopes into the sequence was facilitated using stability and flexibility-inducing linkers, GPGPG. The linkers were instrumental in preserving the epitopes’ structural stability and promoting their pliability. The MHC-II-restricting epitopes were also conjugated with the GPGPG linkers that featured a disulfide cleavage site, thereby enabling the accurate presentation of the epitopes via MHC class I molecules. The MHC-II-specific epitopes were conjugated with linkers with a disulfide cleavage site designated as AAY, facilitating their optimal display by MHC class II molecules. The sequence was finalized by incorporating histidine tags, which have the potential to facilitate the purification and detection procedures of the construct.

Figure 4 The designed vaccine construct and its characteristics. (A) The final vaccine construct; (B) PSIPRED annotation of the secondary structure indicating helices (pink), coils (gray), and sheets (yellow); (C) MEMSTAT prediction of vaccine localization in the cell by amino acids; (D) MEMSTAT depiction of vaccine localization in the cell; (E) disorder prediction by amino acids; (F) The RAMA plot; and (G) 3D vaccine construct refined by Galaxy REFINE.

Immunogenicity and physiochemical assessments

The vaccine construct exhibited significant antigenicity, as evidenced by a score of 0.8742, which suggests its capacity to elicit an immune response. Additionally, it was shown to be non-toxic and nonallergenic, proving its safety profile. The design consisted of 395 amino acids and was anticipated to be electrically neutral at an isoelectric point of 5.41. The vaccine candidate showed an estimated half-life of 1 h in human reticulocytes, indicating stability and predicting quick bodily removal after delivery. Additionally, the vaccine construct’s solubility was shown with a score of 0.637, consistent with its anticipated hydrophilic nature. This attribute proves the construct’s capacity to dissolve and disperse in aqueous settings. Supplementary Table S4 provides further information about the remaining physicochemical properties of the vaccine construct, enhancing the thoroughness of the composition and behavior study.

Homology and human uptake analysis

The human gut microbiome, mouse, and human proteomes, and known human anti-targets showed no substantial sequence similarity to the vaccination design. This data suggests that if administered, it is unlikely that the vaccine cause unexpected effects. By showing a 58.8% sequence similarity to B. licheniformis in the BLAST analysis, the vaccine candidate construct was confirmed to have kept the distinctive qualities necessary for accurately identifying its original proteins and organisms. The specificity and originality of the construct are further highlighted by this observation, improving its potential as a strong vaccination candidate. Additionally, we discovered striking parallels between the vaccine’s design and numerous entries in the Drug bank database.

Secondary structure prediction

Several prediction algorithms were used to analyze the vaccine candidate’s structural properties. According to the PSIPRED analysis, the bulk of protein’s amino acids produced coil structures, with a lesser proportion helping to create beta sheets and alpha helices (Figure 4B). Additionally, MEMSTAT analysis revealed that the vaccine candidate’s majority of residues were projected to be cytoplasmic, indicating the possibility of cellular penetration (Figure 4C). The 338–368 and 372–387 amino acid residues were discovered to span the membrane, whereas 368–372 residues were projected to be extracellular, suggesting the TM nature of the protein, according to the prediction of membrane-spanning regions (Figure 4D). Protein disorder was projected to be quite low in the vaccine candidate by PSIPRED, with just 4% of residues labeled disordered using a threshold value of 0.5 (Figure 4E). The SCRATCH protein prediction predicted the vaccine candidate’s structure disorder as 5.5%.

Immune simulation

Following the injection, the concentration of antigens exhibited a swift increase, reaching approximately 700,000 antigen counts per milliliter. On the initial day, the observed counts exhibited a marked reduction, consistent with the projections derived from the half-life assessment. The data presented in Figure 5A indicate an inverse relationship between the counts of immunoglobulin M (IgM) and immunoglobulin G (IgG) antibodies and antigen counts. Specifically, as antigen counts decreased, the combined counts of IgM and IgG antibodies increased, reaching a maximum of approximately 8500 counts on day 13 of post-injection. Subsequently, there was a gradual decline in the counts of IgM and IgG antibodies, with a decrease to around 3000 counts observed on day 35 post-injection. The concentration of IFN-γ peaked at 420,000 ng/mL and exhibited stability during the period spanning from day 12 to day 15 post-injection. In addition, a limited number of cytokines, including transforming growth factor-beta (TGF-β), interleukin 10 (IL-10), and interleukin 12 (IL-12), were generated. The vaccine administration elicited a minor danger signal, resulting in a surge of IL-2 concentrations, that is, up to 300,000 ng/mL, with the highest point observed on day 5 post-injection, as depicted in Figure 5B.

Figure 5 Immune simulations of the target vaccine candidate. (A) Antigen vs antibody count; (B) chemokine and cytokine counts; (C) B-cell population count; (D) B-cell population per state count; (E) T-helper (TH) cell population count; (F) TH cell population count (in percentage and cell/mm³); (G) natural killer (NK) cell population count; and (H) macrophage (MA) cell population count. Notes: Act=active, Intern=internalized the Ag, Pres II= presenting on MHC II, Dup = in the mitotic cycle, Anergic = anergic, Resting = not active.

A steady rate of 225/mm³ of B memory cells was consistently observed during the production process until day 35 of injection. The generation of non-memory cells exhibited a satisfactory outcome, with an initial decline followed by stabilization at 100 cells/mm³, as illustrated in Figure 5C. Figure 5D demonstrates a rise in active B-cells to 500 cells/mm³ on day 6 post-injection, which remained stable afterward. The data in Figure 5A indicate that the number of presenting cells positively correlated with the increasing levels of antigens. However, on day 8 of the injection, the number of presenting cells decreased to zero (Figure 5D). Significantly, the vaccine construct demonstrated efficacy, as there was insignificant production of anergic cells for the first 5 days only. Figures 5E–5H further depict T-cell populations and other pertinent simulations.

Tertiary structure prediction and validation

Five models were generated using GalaxyRefine to refine our vaccine construct, which 3DPro predicted initially. Model-1 (Figure 4G) was selected from the set of models for validation and subsequent analysis. The preliminary configuration demonstrated a clash score of 13 and a RAMA-favored score of 93.9%. Following refinement, the RAMA-favored value exhibited a notable increase to 94.4%, signifying an enhanced conformational state. Additionally, the clash value decreased to 10.02, indicating a reduction in steric clashes. The RAMA plot findings (Figure 4F) indicate that a significant proportion of amino acids were situated in the highly preferred area, constituting 96.3% of the overall amino acid count. Furthermore, the analysis revealed that 2.3% of amino acids were situated within the additionally permitted region, suggesting a degree of adaptability in those residues.

Epitope validation analyses

As depicted in Figure 6A, the analysis revealed that the vaccine candidate underwent proficient proteasomal cleavage at a critical point of 0.3, evaluated by computerized software. The computational tool ELLIPRO predicted the existence of 11 epitopes, comprising five linear and six discontinuous conformational epitopes. The epitopes were depicted through a visual representation utilizing distinctively colored surfaces. The epitopes identified as B-cell epitopes (Figure 6B) were present in the predicted epitopes, along with other epitopes initially considered for T-cell stimulation. This observation underscores the potential of vaccine candidate to activate both B-cell and T-cell immune responses.

Figure 6 Epitope validation and molecular docking analysis. (A) Vaccine’s immunogenicity prediction by NetChop. Regions in green indicate potential immunogenicity, while regions in pink indicate no immunogenic potential; (B) linear B-cell epitopes predicted in the final vaccine construct. The epitopes are shown in surface of varying colors; docked complexes of: (C) CD4-EPITOPE-MHCII; (D) CD8-EPITOPE-MHCI; (E) TLR1 vaccine; (F) TLR2 vaccine; and (G) TLR3 vaccine.

Molecular docking

In all, 10 putative complexes were generated for each docking pair using the HDOCK server, and the optimal model was subsequently chosen for further analysis. The best models, as depicted in Figure 6, exhibited the least energy score and bond energy. The results presented in Supplementary Table S5 provides the best models based on the highest docking energy. Within these, it is observed that the scores of -244.54 kcal/mol and -231.06 kcal/mol exhibit the most significant molecular interactions between the complexes of MHC class I epitopes and MHC class II epitopes, respectively, with their corresponding epitopes. Figures 6A and 6B illustrate the docking of the optimal complex formed between MHC class I/epitope and MHC class II/epitope with CD8 and CD4 receptors with energy scores of -287.48 kcal/mol and -240.96 kcal/mol, respectively. Figures 6C–6E depict molecular interactions between the vaccine construct and TLRs 1, 2, and 5 with docking energies of -258.49, -239.12, and -275 kcal/mol, respectively.

Molecular dynamics simulations

The MD simulations on all docking complexes revealed that their structures were flexible and less deformable. Figure 7 displays MD simulations of complexes between our vaccine candidate and TLRs. The complexes were stable, as indicated by B-factor and deformability graphs (Figures 7A–F). The eigenvalues calculated for these complexes were 1.690×10^-6 for the TLR-1/vaccine complex, 3.425×10^-7 for the TLR-2/vaccine complex, and 3.612×10^-7 for the TLR-5/vaccine complex. The eigenvalues provided additional verification regarding the stability and adaptability of docked complexes. Covariance maps were produced (as shown in Figures 7G–7I) to detect areas of correlated (illustrated in red) and anti-correlated (illustrated in blue) amino acids within the complexes. The cartographic representations emphasized the areas of the highest correlation among amino acids and yielded valuable perspectives on the structural dynamics of the complex. Furthermore, the elastic network analysis (Figures 7J–7L) demonstrated minimal disruption in atomic correlations, indicating strong and stable molecular interactions within the complexes.

Figure 7 Molecular dynamics simulations of TLR vaccine-docked complexes. (A–C) B-factor graphs of TLR1-; TLR2-, and TLR3-vaccine complexes, respectively; (D–F) deformability graphs of TLR1-; TLR2-, and TLR3-vaccine complexes, respectively; (G–I) covariance maps of TLR1-; TLR2-, and TLR3-vaccine complexes, respectively; (J–L) elastic network maps of TLR1-; TLR2-, and TLR3-vaccine complexes, respectively.

Limitations of computational models and the need for experimental validation

While computational analyses provided valuable insights, they cannot fully capture the complexity of biological systems. Predictions of antigenicity, allergenicity, toxicity, and MHC binding remain dependent on algorithms and limited allele coverage, introducing uncertainty into the results. Therefore, urgent experimental validation is essential. In vitro assays (expression, antigen–antibody reactivity, and cytotoxicity) and in vivo studies in animal models are required to confirm the construct’s immunogenicity, stability, and safety. The absence of an adjuvant, while simplifying design, may also limit immune potency and should be explored in the future experiments. Until such validations are performed, the present findings should be considered preliminary but promising.

Discussion

This study aimed to develop a computational vaccine directed toward B. licheniformis, a pathogen that threatens the well-being of both humans and animals. Infections caused by B. licheniformis are associated with considerable morbidity and economic impact, resulting in various health complications and financial losses.^5,10 Using a comprehensive computational strategy, this study identified potential vaccine candidates and evaluated their interactions with key immune receptors through molecular docking. The 16S rRNA gene was amplified by PCR using specific primers,^19,20 containing conserved and variable sections for species-level identification. This provided genetic confirmation of the isolated strain and guided the downstream vaccine design process.⁴³ Comparative genomic and proteomic approaches subsequently identified unique proteins of B. licheniformis with potential immunogenicity.

Among the screened candidates, proteins of approximately 110 kDa were prioritized, consistent with previous studies,⁴⁵ and only non-homologous proteins with respect to the human proteome, gut microbiome, and known anti-targets were retained to maximize safety. This selective approach ensured that the designed vaccine construct avoided potential off-target effects, a critical factor for therapeutic development.⁴⁶ In addition, proteins lacking extensive TM helices were chosen, as these improve the feasibility of downstream experimental validation.⁴⁷

Docking analyses provided a key finding of this study: the designed vaccine showed strong binding affinities with immune receptors CD4 and CD8 as well as TLRs (TLR1, TLR2, and TLR5), suggesting its ability to trigger both innate and adaptive immune responses.^32,48-49 This finding highlights the construct’s potential to induce a broad and effective immunological response against B. licheniformis. Druggability assessment further reinforced this potential, as the vaccine construct showed significant similarities with multiple entries in the DrugBank database,^50,51 pointing to its translational applicability.

Another important design decision was the omission of an adjuvant. While adjuvants enhance immune responses,^17,52 they often pose formulation and regulatory challenges.⁵³ By excluding adjuvants, the immunogenicity observed in silico can be attributed directly to the construct itself, eliminating confounding factors.⁵⁴ Similarly, the inclusion of linkers ensured structural stability and proper epitope orientation, consistent with evidence that linker-based constructs enhance immune recognition.^55,56 This design choice is supported by earlier findings where linkerless constructs failed to elicit sufficient antibody responses.⁵⁷

Despite the encouraging in silico results, this study has several limitations. The vaccine construct was designed and validated exclusively through computational methods, which, although powerful for early-stage screening, cannot fully replicate biological complexity. The predicted antigenicity, safety profile, and immune interactions therefore require experimental verification. In vitro assays (e.g., protein expression, antigenicity, immunogenicity, and cytotoxicity tests) and in vivo studies in suitable animal models are needed to confirm the construct’s stability, immunogenic potential, and safety. Moreover, the absence of an adjuvant in this design, while simplifying the formulation, may restrict the magnitude of immune responses, which should be addressed in the future wet-lab studies. Addressing these aspects through laboratory experimentation is essential before this computationally designed vaccine can be advanced toward preclinical and clinical development.

In conclusion, this study successfully identified and designed a multi-epitope vaccine construct against B. licheniformis that demonstrated strong in silico interactions with major immune receptors, fulfilled essential safety criteria, and showed promising druggability potential. While experimental validation is still required, these findings provide a solid foundation for the future evaluation and contribute to the broader effort of developing vaccines under the One Health framework, addressing both human and animal health challenges.

Conclusions

This study successfully isolated and identified B. licheniformis, revealing significant similarities with established strains through 16S rRNA analysis. A pan-genome analysis identified 1476 unique proteins, with eight selected as potential vaccine candidates based on epitope prediction, which covered 89.42% of the population. The developed vaccine construct demonstrated substantial antigenicity, stability, and safety, with immune simulations indicating a notable increase in antibody levels. While these findings highlighted the potential of B. licheniformis as a viable candidate for vaccine development, the study’s reliance on computational approaches highlights the need for subsequent validation through in-vitro and in-vivo studies. These validation studies are critical to confirm the immunogenicity, safety, and efficacy of the proposed vaccine in real-world scenarios. This study lays a strong foundation for vaccine development against B. licheniformis by identifying promising candidates and providing preliminary evidence of their potential. However, translating these findings into practical applications will require rigorous experimental research and interdisciplinary collaboration to address the identified limitations and fully realize the vaccine’s translational potential.