Four antigen proteins, namely, A5L, A15L, A35R, and B6R, were chosen according to their function in the life cycle and immunogenicity of the mpox virus. This selection also aimed to assess the immune responses generated by the vaccine against this combination of antigens. They are the orthologs of vaccinia virus glycoproteins—A4L, A14L, A33R, and B5R, respectively (Shchelkunov et al., 2002). These antigens are highly conserved among the mpox, vaccinia, and variola viruses. The protein sequences of mpox virus antigens reported on the National Center for Biotechnology Information (NCBI) database were retrieved. The number of downloaded sequences for A5L, A15L, A35R, and B6R was 1819, 355, 349, and 376, respectively. All the sequences of A5L were used for further analyses as they met the selection criteria such as the known length of the protein sequence. For the rest of the three proteins, the redundant sequences were removed according to their length and quality−179 sequences of A15L, 178 sequences of A35R, and 178 sequences of B6R were used for further analyses. The software Jalview 2.11.1.4 was used for the alignment and conservancy analysis (Waterhouse et al., 2009). The threshold filter of conservation was adjusted to 10 below the threshold, and the consensus sequence was used for B-cell and T-cell (MHC-I and MHC-II binding) epitope prediction.

NetMHCpan EL 4.1 server was applied for the prediction of cytotoxic T lymphocyte (CTL) epitopes (https://services.healthtech.dtu.dk/service.php?NetMHCpan-4.1) (Reynisson et al., 2020). This server predicts the binding of peptides to any MHC-I molecule of the known sequence using artificial neural networks (ANNs) that are a collection of simple interconnected algorithms processing information according to the external input. The globally most prevalent alleles HLA-A^*01:01, HLA-B^*07:02, and HLA-B^*08:01 were selected from the “Allele Frequency Net Database” and used for the prediction of the MHC-I binding epitope. The default parameters of the server were used: the threshold for strong binders was set to 0.5 (% rank), while for weak binders it was 2%. The predicted peptides were identified as strong binders when the % rank was below the set threshold (0.5%), while the weak binders were considered when the peptide rank was above 0.5% and below 2%. For the prediction of helper T lymphocyte (HTL) epitopes, the NetMHCIIpan-4.0 server which performs the prediction of peptide binding to any MHC-II molecule of the known sequence using ANNs was used. The worldwide most prevalent MHC-II alleles DRB1^*07:01 and DRB1^*15:01 were selected for the prediction of MHC-II binding peptides of the mpox virus conserved sequences. According to the default parameters, the predicted peptides were identified as strong binders when the % rank was below the set threshold (1%), while the weak binders were considered when the peptide rank was above 1% and below 2%. The epitopes with the strongest binding capacity to MHC-I and MHC-II alleles were selected. Both NetMHCpan-4.1 and NetMHCIIpan-4.0 use the NNAlign_MA machine learning framework, which has ANNs at its core. NetMHCpan-4.1, a computational tool, has been trained using a comprehensive dataset comprising over 850,000 peptides with quantitative binding affinity (BA) and mass spectrometry-eluted ligand (EL) information. The BA data utilized in the training of NetMHCpan-4.1 encompass 170 MHC molecules from various species—humans (HLA-A, HLA-B, HLA-C, and HLA-E), mice (H-2), cattle (BoLA), primates (Patr, Mamu, and Gogo), equine (Eqca), and swine (SLA). The EL data used in the analysis include 177 MHC molecules from a range of species, such as humans (HLA-A, B, C, E), mice (H-2), cattle (BoLA), primates (Patr, Mamu, and Gogo), swine (SLA), equine (Eqca), and dogs (DLA). NetMHCIIpan-4.0 has been trained using a comprehensive dataset comprising more than 500,000 measurements of BA and EL. This dataset covers the three human MHC class II isotypes, namely HLA-DR, HLA-DQ, and HLA-DP, along with the mouse molecules (H-2). The inclusion of EL data expands the coverage of MHC-II molecules, as the BA data cover 59 molecules, while the EL data cover 74 molecules. Both NetMHCpan-4.1 and NetMHCIIpan-4.0 are the currently recommended algorithms in Immune Epitope Database (IEDB). They demonstrate cutting-edge performance and surpass other algorithms in benchmark experiments (Reynisson et al., 2020).

For predicting the linear B-cell epitopes, the IEDB with the method of Bepipred linear epitope prediction tool (v2.0) was utilized (Jespersen et al., 2017). The BepiPred-2.0 server performs the prediction of B-cell epitopes from a protein sequence via a random forest algorithm trained on epitopes and non-epitope amino acids determined from crystal structures. The threshold was set at 0.6 meaning that the residues scoring above the threshold were predicted to be part of the epitopes. The B-cell epitopes were then selected.

After the prediction of B-cell, MHC-I, and MHC-II epitopes, the epitopes with the most optimal features were selected. Multi-peptide vaccine construct was designed and epitopes of A5L, A15L, A35R, and B6R were fused via the appropriate linkers. The linkers play a crucial role in multi-epitope vaccine development as their reasonable selection allows to enhance structural stability, flexibility, and proper folding, as well as to increase the immune response by allowing to include multiple epitopes in a single vaccine construct. Thus, the well-described linkers KK, GPGPG, and EAAAK were used to connect different peptide components. EAAAK is a rigid linker that provides an alpha helix-forming structure between the domains that enhances the stability and maintains the constant spacing between the domains of the fusion protein as well as preserves their functions (Chen et al., 2013). B-cell epitopes were linked with a flexible KK linker (Rahmani et al., 2021; Tarrahimofrad et al., 2021). The KK linker allows the presentation of each linked epitope to antibodies while avoiding antibody induction against the whole joined sequences resulting in the antibody reactivity to each of the B-cell epitopes (Yano et al., 2005). MHC-I and MHC-II epitopes were conjugated together with another flexible linker GPGPG which is a glycine-rich linker, in addition enhancing the solubility of vaccine construct, and also provides high accessibility and flexibility for the adjacent epitopes (Dong et al., 2020; Khan et al., 2021; Martinelli, 2022). Except for the epitopes, to rationally optimize the fused multi-peptide construct for the optimal immunogenic outcome, adjuvants were also fused. Particularly, 45 amino acid long human β-defensin 3 (hBD3) with immunomodulatory effect (Hoover et al., 2003; Ali et al., 2017; Qamar et al., 2020) was adjoined to the N-terminal of the vaccine after the leading sequence of tissue plasminogen activator which also facilitates the antigen presentation (Ahammad and Lira, 2020). hBD3 is one of the top five adjuvants utilized for COVID-19 subunit vaccines (Mekonnen et al., 2022). Furthermore, defensins play a crucial role in defending against pathogen infections as they effectively bridge the innate and adaptive immune responses through leukocytes such as dendritic cells (DCs) and T cells (Bellamkonda et al., 2022). The rigid linker EAAAK (Nezafat et al., 2014; Tarrahimofrad et al., 2021) was used to conjugate the hBD3 with the pan-HLA-DR-binding epitope (PADRE) which triggers antigen-specific CD4+ T cells (Alexander et al., 2000). EAAAK was used to connect the PADRE sequence to the rest of the construct as well. For future in vitro validation of the expression of the multi-epitope construct, the C-terminal polyhistidine tag-−6xHis—was also adjoined (Boyoglu-Barnum et al., 2021). The untranslated regions (UTRs) were included for ensuring the proper stability and translation of the multi-epitope construct when it is transcribed into mRNA and for allowing the translation efficiency and stability of mRNA while validating the protein expression in vitro and in vivo. The following UTRs were used: cytomegalovirus (HCMV) immediate early (IE) partial sequence as 5′UTR, and human growth hormone partial sequence as a 3′UTR as they have demonstrated the effective stimulation of the protein production in vivo (Rybakova et al., 2019).

The vaccine should have wide human population coverage. Thus, the IEDB population coverage analysis tool was employed to evaluate the worldwide human population coverage of the designed DNA vaccine (Bui et al., 2006). The allergenicity of the designed vaccine was predicted via the AllerCatPro server which predicts the allergenic proteins based on the similarity of both their amino acid sequences and 3D structures compared with the comprehensive dataset of already known allergens from the WHO/International Union of Immunological Societies, Comprehensive Protein Allergen Resource, Food Allergy Research and Resource Program, UniProtKB, and Allergome (Maurer-Stroh et al., 2019; Nguyen et al., 2022). VaxiJen v2.0 was applied to predict the antigenicity of the multi-epitope protein construct (Doytchinova and Flower, 2007a,b). The threshold was set as 0.4, and default parameters were employed (Ghafouri et al., 2022).

The multi-epitope vaccine physico-chemical properties were assessed by ProtParam (Gasteiger et al., 2005). Molecular weight (MW), the composition of atoms and amino acids, instability, estimated half-life, theoretical isoelectric point (pI), aliphatic indexes, and grand average of hydropathicity (GRAVY) of the vaccine construct were evaluated. The MW was determined by summing amino acids' average isotopic masses and the average isotopic mass of a single water molecule. pI was calculated by amino acids' pKa value. Half-life estimation describes the necessary time from protein synthesis until the disappearance of its half-amount from the cell. The instability index was used to determine how stable is the protein in a test tube. The value <40 estimates the protein as stable. GRAVY calculates the hydropathicity index of all amino acids divided by the length of the sequence. The larger number denotes more hydrophobicity of amino acids (Gasteiger et al., 2005). An aliphatic index is used to characterize the relative volume of protein which is occupied by aliphatic side chains and is considered a positive factor in increasing thermostability. The solubility of the vaccine protein was predicted via DeepSoluE. This server is a deep-learning predictor, and its prediction algorithm outperforms the other servers for the prediction of protein solubility (Wang and Zou, 2023).

The 3D structure of the designed vaccine was constructed via RoseTTAFold (Baek et al., 2021). This tool uses a three-track network: integrating protein sequence patterns, amino acid interactions, and 3D structure to accurately predict protein structure and interaction (Baek et al., 2021). After the 3D structure of the vaccine protein was predicted, it was optimized with GalaxyRefine and FG-MD which is a molecular dynamics (MD)-based algorithm for protein structure refinement at the atomic level (Zhang et al., 2011; Heo et al., 2013). This immunoinformatics approach uses molecular dynamics simulation and relaxes and refines the protein structure (Heo et al., 2013). Predicting and knowing the potential structure is important as it determines the function of the protein that allows its application. As the predicted 3D structure provides the understanding of its interactions with the immune system, it is then used for further analyses, e.g., prediction of B-cell discontinuous epitopes and molecular docking. The tertiary structure quality of the 3D protein was then verified with the Ramachandran plot and ERRAT. The Ramachandran plot illustrates energetically permissible regions for the backbone dihedral angles ψ and ϕ of amino acid residues within a protein structure allowing the assessment of 3D protein structure quality (Nelson et al., 2021). ERRAT also represents the quality but for non-bonded interactions, and a higher score of ERRAT indicates a higher quality of the tertiary structure of the protein.

To predict the discontinuous B-cell epitopes, ElliPro was employed. This tool is the recourse of the immune epitope database (IEDB). It allows for predicting the conformational B-cell epitopes according to the antigen's tertiary structure. ElliPro associates predicted discontinuous epitopes with protrusion index (PI) which is averaged over epitopes' residues. The approximation of the antigen's tertiary structure is achieved via the number of ellipsoids. The PI score is defined based on the residue's mass center which is lying outside the largest ellipsoid. A higher value suggests higher solvent accessibility which, on the other hand, is crucial in protein stability and folding (Ponomarenko et al., 2008).

The C-ImmSim server was utilized to analyze the simulated humoral and cellular immune responses elicited against the multi-epitope vaccine protein. After submitting the sequence of vaccines in the FASTA format, the server predicts the immune responses via a position-specific scoring matrix and machine learning algorithms (Rapin et al., 2011). Two HLA-A (01:01 and 02:01), two HLA-B (07:02 and 08:01), and two DRB (07:01 and 15:01) were selected through Allele Frequency Net Database (AFND). The prime, second, and booster doses of the vaccine were injected with 4 weeks gap. The volume and steps of the simulation were 10 and 270, respectively. The vaccine did not contain LPS. The random seed was 12,345, and the injection time periods were 1, 85, and 169 (Rapin et al., 2010).

To evaluate the binding capacity of the vaccine antigen to its recognition receptors of the host, molecular docking analysis was conducted by the server of ClusPro. ClusPro determines the molecular docking of two proteins providing an outcome via presenting putative protein complexes in a ranked list. Ligand conformation, orientation, and position along with the assessment of binding affinity are the major properties that determine molecular docking. Ultimately, the electrostatic-favored protein–protein docked complexes with favorable desolvation-free energies are picked (Kozakov et al., 2017). Analysis of vaccine antigen's docking with MHC-I (HLA-A^*02:01) (6TDS) and MHC-II (HLA-DRB1^*01:01) (1AQD) host receptors was performed. The MHC molecules were used to dock the following epitopes in the vaccine construct: MHC-I ligands—ETNDLVTNVY (A5L), SPTRTWKVL (A15L), DSDVSQEVRKY (A35R), and VSDYVSELY (B6R); MHC-II ligands—APILLPSSTAPVLKP (A5L) and HSDYKSFEDAKANCA (A35R). The tertiary structure of these epitopes was predicted with AlphaFold2 (Jumper et al., 2021; Mirdita et al., 2022). The protein databank (PDB) files 6TDS (MHC-I) and 1AQD (MHC-II) were edited and cleaned to remove heteroatoms, bound peptides, and water molecules. The server of ClusPro performs the docking of rigid body via sampling billions of conformations, clusters 1000 lowest energy structures to find the most favorable models, and refines selected structures via minimization of energy (Kozakov et al., 2017).

All the predicted B-cell, MHC-I, and MHC-II epitopes are given in Supplementary materials 1–12. In total, 30 epitopes were used for the final vaccine construct design. The open reading frame (ORF) of the vaccine construct is shown in Supplementary material. The selected four antigens of the mpox virus are shown in Figure 1. The location, functions, and immune properties of the selected antigens are listed in Table 1. The strategy of in silico design and evaluation of a potentially universal DNA vaccine for mpox are shown in Figure 2. All the epitopes included in the proposed multi-epitope vaccine construct along with the DNA vaccine are given in Figure 3.

The population coverage analyses demonstrated that the 5 MHC-I alleles (HLA-A^*01:01, HLA-A^*02:01, HLA-B^*07:02, HLA-B^*08:01, and HLA-A^*02:12) with IC50 binding affinity <500 nM with the MHC-I epitopes and the 13 MHC-II alleles (HLA-DRB1^*07:01, HLA-DRB1^*01:01, HLA-DRB1^*09:01, HLA-DRB1^*10:01, HLA-DRB1^*04:01, HLA-DRB1^*15:01, HLA-DRB1^*14:01, HLA-DRB1^*13:01, HLA-DRB1^*08:01, HLA-DRB1^*03:01, HLA-DRB1^*11:01, HLA-DRB1^*12:01, and HLA-DRB1^*16:01) with IC50 binding affinity <500 nM with MHC-II epitopes cover 95.21% of the human population, which means that these epitopes could cover approximately 95% of the population from different regions of the world (Supplementary Figure 1 and Supplementary Table 1).

The predicted 3D structure of the vaccine construct along with the quality evaluation is given in Figure 4. The distribution of amino acid residues in the Ramachandran plot was as follows: refined structure had 279 residues in the most favored region, 93 residues in the additionally allowed region, 13 in the generously allowed region, and only 7 residues in the disallowed region. The Ramachandran plot allows us to visualize the energetically allowed and disallowed regions for the dihedral angles. For example, in poor quality homology models, many dihedral angles can be found in the forbidden regions of the plot which indicates the problems with the structure. The favored regions in the Ramachandran plot correspond to the regular secondary structures. The ERRAT score was 88.7218 meaning that all these parameters indicate the favorable overall quality of the vaccine antigen protein.

The instability index of the proposed vaccine construct was 29.12, which classifies the vaccine as stable. The aliphatic index was 60.22 which indicates the high thermostability of the vaccine protein. GRAVY was −0.678, stipulating hydrophilicity. The half-life in vitro was estimated to be 30 h, while the in vivo in the yeast and E. coli were estimated to be over 20 h and over 10 h, respectively. The number of amino acids was 542, MW was 56.55076 kDa, and the theoretical pI was 9.28. The vaccine construct contains 46 negatively charged amino acid residues (Asp+Glu) and 64 positively charged amino acid residues (Arg+Lys). The probability of the vaccine solubility was 0.7963, indicating the solubility of the vaccine protein (Wang and Zou, 2023).

All three conformational epitopes along with the scores are listed in Table 2. The locations of each discontinuous epitope in the tertiary structure of the antigen protein are illustrated in Figure 5.

A vaccine protein construct did not show evidence of allergenicity allowing us to presume that the multi-epitope construct can be considered with low allergenic potential. On the other hand, the multi-peptide was demonstrated to be a probable antigen with a score of 0.4248.

The molecular docking with MHC molecules demonstrated stable binding. ClusPro identifies the most probable complex models by determining the largest clusters available. The interaction positions of the docked complex of vaccine epitopes and MHC molecules are shown in Figure 6. The interaction energy is given in Table 3 where the cluster size is the size of the largest cluster found, while the lowest energy is the Gibbs free energy of the best complex conformation (Kozakov et al., 2017).

In silico analyses demonstrated that the immune response triggered by the proposed vaccine was compatible with the responses generally induced via in vivo immunization. The immune responses after the additional two booster doses of vaccine were stronger compared to the prime immunization. The immunization lowered the antigen level, while high levels of antibodies (IgG and IgM) were produced. Immunoglobulin and antigen levels varied over time. The antigen abundance peaked at each injection time point (Figure 7A). The humoral response was induced after each shot, and the antibody levels remained elevated during the weeks after the last vaccine injection. The immunization stimulated the B cells (Figures 7B, C). The ignorable number of anergic cells and activation of CTLs was observed (Figure 7D). The counts of HTLs demonstrate that the duplication phase starts immediately after each injection (Figure 7E). Days after each injection and upon successful interaction with antigen-presenting cells, HTLs start to duplicate and differentiate into memory cells. HTLs also foster the secretion of cytokines and humoral response. High levels of IFN-γ, TGF-β, and IL-10 were induced upon each injection (Figure 7F). The other parameters of immune responses are shown in Supplementary Figure 2.