A comprehensive collection of 198 MAP isolate genomes was sourced from NCBI. The CDSs of all genomes underwent annotation with Prokka (version 1.14) (37). Utilizing Roary (version 3.13.0) (38), the proteomes of 198 samples were analyzed with the aim of determining the core proteomes; this was done using the default settings.

Following the identification of core proteomes through pan-genome analysis, subcellular localization was predicted using the PSORTb server (39). Subsequently, the antigenicity of the curated protein amino acid sequences was meticulously evaluated using VaxiJen (40).

To ensure the efficacy and safety of our multi-epitope vaccine, we leveraged various bioinformatic tools to thoroughly evaluate and refine our epitope selections. CTL epitopes were predicted using the NetCTL 1.2 server, which supports predictions for CTL epitopes restricted to 12 MHC class I supertypes: A1, A2, A3, A24, A26, B7, B8, B27, B39, B44, B58, and B62. MHC class I binding and proteasomal cleavage employ artificial neural networks, while TAP transport efficiency is forecasted using a weight matrix (41). Peptides were predicted to bind to MHC class II molecules using the IEDB MHC class II binding prediction tool, based on the NetMHCIIpan-4.1 method and including HLA class II alleles from three supertypes: HLA-DR, HLA-DQ, and HLA-DP. B-cell epitopes were identified using the ABCpred server (41–43). All predicted epitopes were then evaluated for allergenicity using AllergenFP and AllerTOP (44, 45), toxicity using ToxinPred, and antigenicity using VaxiJen (v2.0). This multi-layered screening process ensured that only non-toxic, non-allergenic, and antigenic epitopes were selected for inclusion in our vaccine design.

In the design of the multi-epitope vaccine, we integrated previously identified T-cell and B-cell epitopes. To facilitate efficient separation of these epitopes in vivo and circumvent the likelihood of unintended junctional epitope creation, tailored linkers were employed: AAY for CTL epitopes, GPGPG for HTL epitopes, and KK specifically for B-cell epitopes (46). In addition, the LTB subunit of the heat-stable enterotoxin from UniProt (Code: P0CK94) was used as an adjuvant. LTB is a potent mucosal adjuvant that induces both Th1 and Th2 immune responses, adaptable to different administration routes for tailored immunity. Its pentameric structure facilitates antigen delivery and enhances immune responses at both systemic and mucosal sites (47–50). This adjuvant was conjugated to the N-terminus of the vaccine via the EAAAK linker, Natural linkers like EAAAK are rigid and stable due to their tightly packed backbone and intrasegment hydrogen bonds, which promote α-helical structures, allowing them to act as rigid spacers between protein domains (51).

The designed vaccine underwent further allergenicity evaluation using AllergenFP and AllerTOP (44, 45). The amino acid sequence composition of the multi-epitope vaccine underwent thorough antigenicity analysis utilizing the VaxiJen tool. Additionally, to gain a comprehensive understanding of its physicochemical attributes, the vaccine was evaluated with ProtParam. This analysis imparted crucial insights into the vaccine’s characteristics, including its stability, solubility, and potential immunogenicity (52).

To delve deeper into the secondary structural components of our vaccine construct, we utilized the freely accessible PSIPRED server, an online protein structure analysis tool. By submitting the primary amino acid sequence of the vaccine peptide as input, PSIPRED predicted its secondary structure, offering invaluable insights into the folding patterns and conformational landscapes that underpin the vaccine’s structural integrity and potentially impact its functional performance (53).

To guarantee heightened accuracy in forecasting the protein’s configuration, the initial three-dimensional blueprint, crafted by the Phyre 2 platform, underwent a meticulous optimization process facilitated by the GalaxyRefine system (54, 55). Following this refinement, the PDB documentation embodying the refined structure was subjected to a rigorous assessment of its tertiary structure’s integrity utilizing the SWISS-MODEL server as a tool (56). Additionally, the validity of the protein’s architecture was confirmed through the utilization of the ProSA web-based service (57).

ElliPro (http://iedb.org/ellipro/) was used to identify B-cell epitopes within the protein’s refined 3D structure (58). It assigns a score, the Average Protrusion Index (PI), to each predicted epitope. Notably, the PI threshold for an ellipsoid that encompasses 90% of the protein residues is 0.9, while the remaining 10% lies outside this boundary. Determining this PI threshold relies on an assessment of the deviation from the center of mass for residues that lie outside the largest reasonable ellipsoid. In comparative evaluations, ElliPro has demonstrated superior accuracy in structure-based epitope prediction, evidenced by its Area Under the Curve (AUC)value of 0.732. This high AUC value underscores ElliPro’s reliability andprecision as a tool for predicting protein epitopes.

To mimic the intricate biological milieu encountered by proteins within living organisms, a molecular dynamics (MD) simulation was executed employing GROMACS (59). Furthermore, the durability and physical characteristics of the docked complex were meticulously evaluated through the iMODS web server (accessible at http://imods.Chaconlab.org) using normal mode analysis (NMA). The iMODS analysis divulged pivotal insights into the complex’s structural dynamics, offering a comprehensive assessment encompassing deformability, B-factor analysis, eigenvalue examination, variance measurements, covariance mapping and elastic network modeling.

The designed vaccine was subjected to various docking simulations to evaluate its binding affinity towards human Toll-like receptors, including TLR4. The crystallized structures of TLR4, designated by the Protein Data Bank with the unique PDB code 3FXI. Molecular docking was conducted utilizing the web-based ClusPro tool (60). The ensuing three-dimensional representations of both the Toll-like receptors and the vaccine-TLR4 complexes were rendered visible using PyMOL software. The specific amino acid sequences that mediated the interplay between the tailored vaccine and TLR4 were pinpointed with the assistance of PDBsum, facilitating a comprehensive understanding of their interaction patterns (61).

To evaluate the immunological response elicited by the vaccine, we used C-ImmSim, an agent-based simulator that employs the Celada-Seiden model to simulate the mammalian immune response to the designed vaccine (62). This simulation generates profiles for both humoral and cellular immune responses. The default settings were retained, with modifications made to the “Simulation Steps” set at 600. Additionally, the timing of injections was adjusted, with the first injection administered on day 1, the second on day 30 and the third on day 60.

The Escherichia coli K-12 strain was selected as the host for the in silico cloning process to create the vaccine construct. The codon sequence of the vaccine was optimized using the JCAT codon adaptation tool (63). Subsequently, the sequence was cloned into the pET28a(+) vector, leveraging the capabilities of SnapGene software for seamless insertion (64).

The annotation of 198 complete MAP genomes was displayed using Prokka, followed by the conversion of genomes into GFF3 format files. Subsequently, a pan-genome analysis was conducted utilizing Roary, which revealed a total of 11,930 genes among the 198 genomes. Among the identified genes, 2,312 were categorized as essential core genes. Our comprehensive analysis generated a whole-genome phylogenetic tree, along with a matrix that clearly indicates the distribution of both core and accessory genes across samples. This information is graphically shown in Figure 2A .

A subcellular localization screening, utilizing PSORTb, classified 991 proteins as cytoplasmic, 417 proteins as cytoplasmic membrane-bound, 15 proteins as extracellular, 7 proteins as cell wall-associated, and 395 proteins had an unknown localization. Proteins with extracellular secretion and those in the cell wall were chosen for detailed analysis ( Figure 2B ).

The proteins identified in the preliminary stage underwent scoring by VaxiJen. From among these proteins, 22 were recognized as potential antigens, satisfying the default threshold of 0.4. Upon further screening, 5 proteins ( Table 1 ) emerged with scores exceeding 0.7. These proteins were identified as high-priority candidates and were selected for epitope extraction to construct the multi-epitope MAP vaccine.

Employing the updated NetCTL-1.2 platform, we forecasted CTL epitopes from the selected protein repertoire. Out of the 227 predicted epitopes, a subset of 9 were meticulously chosen, guided by their high binding affinity ratings with MHC I supertypes, alongside comprehensive evaluations of toxicity, allergen potential, and antigenicity. These non-toxic, non-allergenic, and antigenic epitopes, shown in Table 2 , were deemed suitable for inclusion in the vaccine construct.

Using the IEDB server, we predicted HTL epitopes for three HLA supertypes: HLA-DR, HLA-DQ, and HLA-DP. We screened 10 non-toxic, non-allergenic, and antigenic epitopes for inclusion in the vaccine construct, as shown in Table 3 .

To predict B-cell epitopes, the ABCpred server was utilized with an epitope identification threshold of 0.5. Among the 102 epitopes screened based on their amino acid sequence positions, 10 non-toxic, non-allergenic, and antigenic epitopes were selected for inclusion in the final vaccine construct., as shown in Table 4 .

Using established criteria, the vaccine construct was carefully designed by integrating selected epitopes. A total of 10 B-cell epitopes, 10 HTL epitopes, and 9 CTL epitopes were screened. The B-cell epitopes were linked using the KK linker, HTL epitopes with GPGPG, and CTL epitopes with AAY. The LTB adjuvant was conjugated to the N-terminus of the vaccine via the EAAAK linker. The final vaccine construct, consisting of 594 amino acids. A graphical representation of the vaccine is shown in Figure 3 .

The vaccine’s potential for allergenicity was assessed by analysis on both the AllergenFP and AllerTOP servers. Both platforms indicated that the vaccine sequence is likely a nonallergenic. To further confirm its antigenicity, the VaxiJen server was used for reassessment, which produced a score of 1.1949.

The web server ProtPara was utilized to evaluate various physicochemical properties of our designed vaccine. Comprising 594 amino acids, includes 33 negatively charged and 56 positively charged residues, with a molecular weight of 58 kDa. The theoretical PI of the vaccine was predicted to be 9.52, containing 8205 atoms, and the chemical representation was C₂₅₅₁H₄₀₆₆N₇₃₂O₈₄₁S₁₅. The predicted half-life of the vaccine is 30 hours in vitro, more than 20 hours in yeast (in vivo), and over 10 hours in Escherichia coli (in vivo). The instability index of the protein was 27.26, which is considered stable. Notably, its aliphatic index of 58.69 suggests relatively higher thermostability. Additionally, the vaccine demonstrates a hydrophilic nature, with a grand average of hydropathicity (GRAVY) value of -0.341.

The PSIPRED server was employed to predict the secondary structure of the chimeric peptide, which consists of 594 amino acid residues. The prediction revealed a composition comprising 14.14% alpha helix, 23.23% beta-sheet, and a noteworthy 62.63% of coils, shown in Figure 4A . This breakdown hints at a preponderance of coils in the vaccine’s structure, with smaller proportions of alpha helices and beta-sheets. The arrangement of these secondary structures is crucial for determining the vaccine’s stability, solubility, and immunogenic potential.

To refine the three-dimensional (3D) structure of the designed vaccine, the construction process first took place on the Phyre 2 platform, Visualization of·initial vaccine model was achieved through PyMOL, as shown in Figure 4B . An optimal model typically boasts a higher concentration of residues located in the Ramachandran-favored region, along with a lower proportion in the outlier and rotamer regions. A detailed analysis of the Ramachandran plot for the initial tertiary structure revealed that 78.38% of residues fell within the favored region, 7.43% were outliers, and 0% belonged to the rotamer region ( Figure 4C ). The MolProbity and Clash scores for this structure were 2.85 and 41.1, respectively. Additionally, the ProSA tool was used to evaluate the quality of the initial structure, yielding a Z-score of -3.51, as shown in Figure 4D .

The vaccine’s initial three-dimensional structure was refined using the GalaxyRefine server. Among the 5 generated models, the one boasting a Rama favored score of 81.1 was chosen for deeper analysis. Visualization of·the selected model was achieved through PyMOL, as shown in Figure 4E , the vaccine’s finalized structure exhibited improved metrics, with 80.24% of residues in the Ramachandran-favored region, 5.57% as outliers, and 0.95% in the rotamer region ( Figure 4F ). The MolProbity and Clash scores improved to 1.96 and 4.50, respectively. The ProSA Z-score also improved to -3.83, further confirming that the refined structure aligns even more closely with experimentally determined structures ( Figure 4G ).

The ElliPro tool was used to predict the conformational B-cell epitopes in the multi-epitope vaccine. Epitopes with a score of 0.7 or above were chosen for further evaluation. This selection revealed five discontinuous B-cell epitopes, comprising approximately 270 residues with scores ranging from 0.722 to 0.793. The number of residues within each conformational epitope varied between 20 and 52 ( Figure 5 ).

Molecular docking simulations between the designed vaccine and TLR4 were performed using the ClusPro server. Each simulation produced 30 clusters, and the cluster with the most favorable energy score of -1037.4 was selected as the best docking arrangement. To visualize these docking outcomes, we employed PyMOL ( Figure 6A ).

To achieve a deeper comprehension of the vaccine’s interplay with TLR, we analyzed the participating residues and the binding strength of the receptor-ligand interaction through the PDBsum tool. In the vaccine-TLR4 docking complex ( Figure 6B ), the vaccine presented 33 interface residues, while TLR4 contributed 31, with respective interface areas of 1612 Å2 and 1562 Å2. This complex exhibited 2 salt bridges, 18 hydrogen bonds and 183 non-bonded contacts, providing strong evidence of the strong affinity between the vaccine and TLR4.

Molecular dynamics simulations were performed using GROMACS to assess the biological stability of the vaccine and the vaccine-TLR4 complex. For both simulations, the OPLS-AA force field was employed, with the vaccine centrally positioned in a water-filled cubic box and chloride ions added to maintain neutrality.

In the vaccine protein simulation, energy minimization reduced the system’s energy to approximately -14,000,000 KJ/mol ( Figure 7A ), followed by equilibration at 300 K and 1 atm using NVT and NPT ensembles ( Figures 7B, C ). RMSD analysis revealed a stable conformation with an RMSD of about 0.6 nm after 50 ns ( Figure 7D ). RMSF results indicated high stability, with an average RMSF of 0.45 ± 0.3 Å and flexible residues fluctuating around 0.8 Å ( Figure 7E ). The Rg analysis showed a compact and stable structure, with an average Rg of 2.45 ± 0.05 Å ( Figure 7F ). These findings suggest that the vaccine protein is stable and suitable for further experimental validation.

In the vaccine-TLR4 complex simulation, energy minimization lowered the system’s energy to around -9,000,000 KJ/mol ( Figure 8A ). The system was similarly equilibrated at 300 K and 1 atm using NVT and NPT ensembles ( Figures 8B, C ). RMSD analysis showed a consistent conformation with an RMSD of approximately 0.5 nm after 15 ns ( Figure 8D ). The RMSF results averaged 0.5 ± 0.3 Å, with flexible residues fluctuating around 0.8 Å ( Figure 8E ). The Rg analysis indicated a compact and stable structure, with an average Rg of 5 ± 0.1 Å ( Figure 8F ). Overall, the vaccine-TLR4 complex maintained a consistent and stable structure, indicating high stability throughout the simulation.

Normal mode analysis (NMA) was conducted using the iMOD server to evaluate the stability of the vaccine within the vaccine-TLR4 complex. Regions containing hinge points exhibited considerable deformability, as depicted in Figure 9A . The B-factor values showed a direct correlation with the root mean square in the NMA ( Figures 9B, C ). The eigenvalues for the vaccine-TLR4 complex were calculated to be 1.665426 × 10⁵, indicating the energy required for structural deformation, where lower values suggest easier deformation ( Figure 9D ). The relationships between pairs of residues is visualized in a covariance matrix ( Figure 9E ), while the elastic network model depicted interatomic connectivity through springs ( Figure 9F ). These results suggest that the vaccine maintains a continuous interactions with TLR4.

The analysis of simulation results revealed several noteworthy findings. Firstly, there was a significant increase in antibody titers following subsequent immunizations, indicating a strong humoral immune response ( Figure 10A ). Secondly, the secondary immune response showed marked augmentation compared to the primary response, evident in the heightened activity of B cells ( Figures 10B, C ), helper T cells ( Figures 10D, E ), and cytotoxic T cells ( Figure 10F ). Additionally, the rise in dendritic cells ( Figure 10G ) and macrophage ( Figure 10H ) concentrations pointed to efficient antigen presentation by these APCs. Finally, cytokine and interleukin levels ( Figure 10I ) saw a notable surge following the second inoculation. In summary, these results suggest that the vaccine can stimulate high levels of antibodies, activated immune cells, cytokines, and APCs, underscoring its immunogenicity and potential efficacy against the targeted pathogen.

The JCat server facilitated the codon optimization of the vaccine’s DNA sequence, resulting in a CAI value of 1.0 and a GC content of 53.64%, thus ensuring high expression efficiency. With the aid of SnapGene software, the vaccine sequence’s successful insertion between the BamH I and Hind III restriction sites was visualized, appearing highlighted in gray within the constructed vector ( Figure 11 ).