The study starts by selecting aerolysin as the potential candidate to synthesize a MEV against it. The virulence of aerolysin was confirmed via an online tool called the Virulence Factor Database (VFDB) (22). This database consists of thousands of virulence factors associated with different known pathogens. Furthermore, clinical case studies and the role of aerolysin in causing different infections in humans were confirmed by PubMed and UniProt (23). The amino acid sequence of aerolysin was obtained from the UniProt database, and the structure was obtained from the Protein Database (PDB) with ID (1PRE). Obtained sequences were processed for homology checks against Homo sapiens, Lactobacillus rhamnosus, Lactobacillus casei, and Lactobacillus johnsonni by using BLASTp from the NCBI server (24). This was done to check that the aerolysin sequence does not have any similarity with the human and normal flora proteome. If similarity was shown, the study would not proceed because the obtained epitopes may provoke autoimmunity (25). Additionally, a conservancy analysis was conducted on aerolysin using the Blastp tool (24) to assess the occurrence of this virulence factor in other strains of the Aeromonas family. The findings from this analysis will enhance the significance of aerolysin-based MEV if the majority of strains harbor this gene in their genome. Consequently, such a MEV could confer broad immunity against all strains of Aeromonas that possess aerolysin.

The obtained sequence of aerolysin was subjected to epitope prediction. The Immune Epitope Prediction (IEDB) database is a unique platform that consists of millions of known and predicted epitopes related to critical immune cell receptors, i.e., TLR4, BCR, MHC-II, MHC-I, etc. (26, 27). The system operates on machine learning algorithms to predict epitopes from the query sequence, leveraging the chosen alleles as a basis. This methodology allows for the accurate identification of epitopes tailored to specific alleles, enhancing the precision and effectiveness of the MEV candidate. The MEV designed from the predicted epitopes has demonstrated success across multiple vaccine design projects when tested experimentally (28–30). This validation underscores the efficacy and versatility of the MEV approach, further supporting its potential as a promising solution in vaccine development (31). Our study focuses on B-cell epitope prediction, and the obtained epitopes were separately processed for MHC-I and MHC-II epitope prediction. For the prediction of MHC-II and MHC-I epitopes, all known alleles representing the entire human population were selected. This approach ensures that the predicted epitopes encompass a broad spectrum of alleles related to the immune system, thereby enhancing the versatility and applicability of the MEV for experimental approaches (32). The resulting epitopes were selected on the basis of low percentile ranking because a lower ranking score resulted in the maximum binding of the epitopes to immune receptors (33, 34).

The obtained epitopes were subjected to a series of analyses so that filtered epitopes could be obtained and used as part of MEV. Antigenicity and toxicity were checked for each epitope using Vexigen 2.0 (35) and the ToxinPred tool (36). Epitopes with no toxicity and an antigenicity value of less than or equal to 0.4 underwent additional solubility analysis. Each epitope’s solubility was assessed using the Innovagen Peptide Calculator, and those with good water solubility were processed further (37). To check for allergic reactions, filtered epitopes were examined using the Allertop 2.0 tool. Only those epitopes that were negative for allergenicity were further processed (38).

The vaccine construct was assembled by sequentially linking the filtered epitopes. An adjuvant was initially conjugated to the first epitope via an EAAAK linker. Subsequently, GPGPG linkers were utilized to connect each subsequent epitope in a tandem array (39). The purpose of adding an adjuvant was to enable the MEV to elicit a strong immunological response, and the objectives of adding linkers were to stabilize the vaccine and prevent self-complementary binding of its sequences (40). The resulting construct’s physiochemical properties were also anticipated to verify its theoretical PI, solubility in water, stability, and other properties. Immune simulation tests were conducted for the vaccine construct using the c-IMMSIMM tool in order to evaluate immune responses (cellular and humoral) against the vaccine candidate (41).

A three-dimensional structure was predicted for the MEV via an online structure modeling tool known as iTESSOR (42). iTESSOR uses a hierarchical approach that combines threading, structural refinement, and template-based fragment assembly (43). It uses particular methods to further enhance the structure after making a forecast. Because of this combined method, iTESSOR is now at the forefront of protein structure prediction, allowing for reliable structure prediction for a wide range of protein folds and sequences (44). The obtained model was further processed for refinement via the GALAXY refine tool in order to minimize steric clashes among the amino acids and increase the stability of the MEV by transforming coils into suitable helix structures (45). The final stability of the MEV was checked via an online tool known as ERRAT (46) and PDBsum, where a Ramachandran plot was obtained for the predicted MEV (47, 48).

Molecular docking was used to evaluate the MEV’s ability to bind to immune cell receptors. Using the ClusPro server, the vaccine construct was docked with MHC-I (PDB ID: 1I1Y), MHC-II (PDB ID: 1KG0), and TLR4 (PDB ID: 4G8A) receptors (49). The complexes were further visualized via the UCSF Chimera tool (50), and the PDBsum tool was used to analyze the kinds and quantities of interactions among the docked complexes (51). Furthermore, the docked complexes were processed through molecular dynamic simulation (MDS), which is a method employed to comprehend molecular behavior and characteristics at the atomic level. It proves effective in assessing the binding stability of a complex, such as a ligand protein, within a dynamic environment. The AMBER 20 package was used, in which the system underwent 500 steps of steepest descent and 500 stages of conjugate gradient minimization (52). Position constraints with a force constant of kcal mol⁻¹ Å⁻² were used to keep the protein stable. Subsequently, there were 2,000 more steepest descent stages and 2,000 more conjugate gradient minimization phases, respectively. Secondly, the system was heated to the target temperature of 300 K for 20 ps under weak positional constraints on the protein atoms (force constant of 10 kcal mol⁻¹ Å⁻²) and constant volume periodic boundary conditions (NVT). After that, the system was calibrated by running a production simulation for 100 ns while maintaining the same pressure and temperature (NPT) for roughly 40 ns. Using isotropic position scaling and a relaxation duration of 2 ps, an average pressure of 1 atm was maintained. Langevin dynamics with a collision frequency of one ps-1 were used to control the temperature (53). The Particle Mesh Ewald (PME) method was utilized to handle nonbonded interactions and long-range electrostatic interactions, with a 10-Å limit (54). The numerical integration time step was 2 fs, and the SHAKE method was utilized to limit all bonds containing hydrogen (55). VMD and the PTRAJ program from the Amber11 package were used to analyze the simulation results (56). For receptor–binder complex systems, the binding free energy (ΔGbinding) was calculated using molecular mechanics with a generalized Born and surface area solvation (MM/GBSA) approach (57). Throughout the simulation trajectory, 1,000 pictures were captured at 20 ns intervals in order to calculate the MM/GBSA free energy difference.

The designed MEV showing all the important properties of a common vaccine was subjected to in silico cloning (58). This was achieved by performing codon optimization of the MEV via an online tool called Jcat (59). Jcat tool sets the sequence of the MEV according to the codon usage of the expression system used in in silico cloning (Escherichia coli K-12 was used as an expression system). The processed sequence was inserted into a well-known vector named pet28+(a) via the SnapGene tool (60). Due to its T7 promoter, several cloning sites, His-tag fusion (polyhistidine that makes it easier to purify the produced protein), and a selectable marker, this vector is primarily utilized in cloning techniques (61). This process will yield accurate data that experimentalists may use to allow manufacturing at the industrial level.

The VFDB showed that aerolysin is among the most essential virulence factors of A. hydrophila and plays important roles in human-related infections. The sequence was obtained from the UniProt database, and the structural ID for aerolysin was obtained from PDB, as shown in Table 1 .

Each virulence factor was separately aligned via BLASTp tool against the human proteome and important normal flora bacteria, i.e., Lactobacillus rhamnosus, Lactobacillus casei, and Lactobacillus johnsonni, in order to check if there is sequence similarity or not. The sequence of each organism showed no such similarity with the virulence factor, and it was safe to use for further analysis. This was done because, while constructing such a vaccine candidate, one should keep in mind that the sequence used in constructing MEV should not resemble the sequence of these organisms because it may lead to autoimmune disorders or abnormal killing of the normal flora of the human body because our adaptive immunity can be provoked against our own body and normal flora. Secondly, sequence alignment of aerolysin with all known strains of the Aeromonas family reveals that the selected virulence factor is present in all those strains of Aeromonas known for causing serious infections in the human population. Supplementary Table S1 shows that aerolysin is present in 17 strains of A. hydrophila, 39 strains of Aeromonas spp., three strains of A. salmonicida, and each strain of A. bestiarum, A. piscicola, A. caviae, and A. dhakensis, respectively. Alignment scores of the aerolysin with the mentioned strains showed a maximum result of 95%–100%, which clearly illustrates that these genes are exactly present in these mentioned strains.

The chosen proteins were ranked in order of priority for immunological epitope prediction. This was accomplished by first predicting the B-cell epitope and then the T-cell epitope. To find T-cell epitopes, additional processing was done on the B-cell epitopes. B cells, macrophages, and cytotoxic T lymphocytes are all stimulated by helper T lymphocytes. On the other hand, antigens can be directly recognized by cytotoxic T cells (62). Conversely, B cells can develop into plasma cells, which produce antibodies (63). Aerolysin was subjected to B-cell epitope prediction. The sequence was filtered and spliced into 14 linearly predicted B-cell epitopes shown in Table 2 . All the epitopes were also presented schematically in Figure 2 , which showed that epitopes present in the yellow region above the threshold value of 0.5 are considered B-cell epitopes, while sequences below the threshold present in the green region do not belong to the B-cell epitope category (26). Out of 493 amino acids of aerolysin, 262 amino acids (selected from each chain, i.e., A and B of aerolysin) were further processed separately for MHC-II and MHC-I epitope prediction. All 262 epitopes of B cells were used for T-cell epitope prediction. The first MHC-II epitopes were predicted, followed by MHC-I epitope prediction. For MHC-II epitope prediction, the IEDB-recommended method was used for epitope prediction (64), and all the known alleles of MHC-II molecules were selected. This was done because the selection of more alleles results in diverse epitopes that could be present in a high percentage of people. Epitopes with a percentile score of ≤ 10.0 were selected, and the remaining were rejected because a lower percentile score results in more efficient binding. For MHC-I epitope prediction, the same IEDB-recommended method was used (64), and all the known alleles for MHC-I were selected. Two filtration thresholds were set for choosing the best epitopes, i.e., top 1% epitopes and epitopes having a percentile score of ≤ 1.0 were selected. Top 1% was selected because this would filter those epitopes that are common in all the selected alleles, and the lowest percentile score results in better binding with the immune receptors. The obtained epitopes of MHC-II and MHC-I are shown in Table 3 .

Out of a multitude of possibilities, only the most promising epitopes made the cut for our vaccine design. After undergoing rigorous assessments for toxicity, solubility, allergenicity, and antigenicity, only a select few emerged victorious. Table 4 showcases these epitopes, ready to serve as the building blocks of our MEV.

In total, 30 unique epitopes were chosen from the list of combined epitopes after the aforementioned analyses. One of the main problems was resolved by creating a multi-epitope-based vaccination construct by joining different kinds of designated epitopes with certain GPGPG linkers. Furthermore, the EAAAK linker was used to link the epitope peptide and the adjuvant for the cholera toxin B component. Because GPGPG linkers can effectively block junctional folding and initiate an immunological response involving T-helper cells, they were placed between epitopes (65). EAAAK is a stiff, stable α-helical peptide linker with an intramolecular hydrogen bond and a closed-packed backbone. Consequently, in a fusion protein, the EAAAK linker serves as a domain spacer (66). Additionally, linkers facilitate the union of epitopes to form a significant structure with a polytope shape (67). Cholera toxin B was used as an adjuvant because it significantly increased the synthesis of IgA in the mucosa and other immune responses (68). The reason behind this is that it is nontoxic and has the ability to attach itself to the monosialotetrahexosylganglioside (GM1) receptor. This receptor can be found in the cytosols and on the membranes of various cells, such as B cells, macrophages, dendritic cells, gut epithelial cells, and antigen-presenting cells (69). In a similar vein, the adjuvant employed is risk-free and produces strong immune responses that are particular to the antigen with which it is coupled (70). Figures 3A, B and Table 5 present the MEV construct and its significant attributes, respectively.

The immunological reactions to the MEV were analyzed using the C-ImmSim server to ascertain if the epitopes would be sufficient to produce immunity (41). By employing this technique, it is also possible to determine the emergence of immunological interactions between the epitopes and specific targets. The ability of the MEV construct to elicit potent cellular and humoral immune responses is shown in Figure 3C . An increase in the development of adaptive responses, such as IgG and IgM antibodies, was seen in a C-immune simulation analysis carried out 35 days after the human immune system was virtually exposed to the highest dosage of vaccine antigen. Similarly, robust cellular immune responses were evident by significant production of interferon-gamma, interleukin (IL)-10, and IL-2 observed within 5 days postadministration. These findings highlight the MEV’s ability to stimulate potent immune responses, indicating its potential as an effective vaccine candidate (71–73).

The MEV construct was modeled for structure prediction. In order to do this, the final MEV construct’s amino acid was uploaded in FASTA format to the iTESSOR program, and an ab initio modeling technique was used to predict the structure (42). Five different structures were forecast in total, and the one with the highest confidence score was chosen. The modeled structure was refined by reducing steric clashes among the residues, making it stable and suitable for further analysis (46, 74, 75). The stability of the MEV structure was predicted, and it showed acceptable values of stability. Two algorithms were used to predict the stability, i.e., the Ramachandran plot and the ERRAT plot. The Ramachandran plot illustrates the division and placement of the amino acids in the MEV construct into distinct areas, each of which denotes a distinct stability level (76). The four zones in this plot are the most favored areas (red), the additional permitted areas (brown), the generously allowed regions (yellow), and the disallowed regions (pale). Following the placement of the MEV construct residues in these areas, the combined stability was computed. Based on each residue’s phi and psi angles, which are shown on the plot’s x-axis and y-axis, residues are arranged in these regions. The majority of the MEV construct’s residues were found to be in the allowed region (75.5%), which was followed by the generously allowed zone (1.3%), the additional allowed regions (19.8%), and the disallowed region (2.4%) ( Table 6 ). Our MEV construct is stable, as evidenced by the presence of more residues in the permitted regions and fewer in the prohibited zone ( Figure 3D ). The ERRAT score for structure stability also lies in the acceptable range of 70%, respectively.

Docking is an essential method for assessing how well two molecules bind together. The MEV candidate was bound to the most important receptors of the human immune system, i.e., MHC-I, MHC-II, and TLR4 in order to predict its binding efficiency. The pdb structures of the immune receptors were retrieved from the Protein Data Bank (PDB) (77) and were separately docked with the MEV, and binding free energies were calculated where the obtained values clearly demonstrate maximum binding ( Table 7 ). Furthermore, the PDBsum tool was used to verify the kinds and quantities of interactions. The MEV exhibited effective binding with TLR4, MHC-I, and MHC-II, with the highest number of contacts. In addition, hydrogen bonds were observed among the complexes, providing additional evidence of maximal binds. Figure 4 displays the structure of docked complexes as well as the quantity, kind, and number of interacting residues.

The predicted MEV was ultimately subjected to in silico cloning, as it was predicted to be the most promising candidate for evoking immune responses and preventing A. hydrophila infection in individuals. To facilitate cloning, the vaccine construct was optimized using a codon adaptation tool tailored for the E. coli K-12 expression system (59). Results indicated that approximately 95% of the codons were successfully modified to align with the expression system’s codon usage, suggesting that the MEV would be efficiently expressed in E. coli ( Figure 5A ). Codon optimization is crucial because the expression efficiency of codons varies among organisms, depending on their specific codon usage patterns (78). The optimized vaccine construct was then successfully cloned into a specialized vector, pet28a(+), with modified restriction sites as depicted in Figure 5B , respectively (79).

All docked complexes subjected to simulations yielded acceptable results. The root mean square deviation (RMSD) plot in Figure 6A showed linear deviations for each complex, indicating stability upon ligand binding and an absence of extreme deviations across the timeframe. This is a kind of ideal situation for our docked complexes. Root mean square fluctuation (RMSF) results in Figure 6B further clarify that amino acids within the active sites exhibit effective fluctuations conducive to efficient engagement and binding with their respective ligands. Molecular Mechanics Generalized Born Surface Area (MM-GBSA) and Molecular Mechanics Poisson–Boltzmann Surface Area (MM-PBSA) analyses demonstrated the stability of each complex upon binding ( Table 8 ). Notably, the obtained energy values were significantly low, signifying maximal stability and efficient binding potentials of the MEV with the immune receptors.