The amino acid sequences of A35R, B6R and H3L of three strains (Zaire-96-I-16, 2001, GenBank: GCA_000857045.1; MPXV-UK_P2, 2018, GenBank: MT903344.1; 2022 reference strain, GenBank: GCA_014621545.1) were obtained in FASTA format from National Center for Biotechnology Information (NCBI) (https://www.ncbi.nlm.nih.gov/) database. To identify the common and specific sites of the single protein from different strains, we performed alignment of the target protein sequences from three objective MPXV strains using the Clustal Omega program (https://www.ebi.ac.uk/Tools/msa/clustalo/) available on the EMBL-EBI server. The visualization of the results of sequence alignment performed by Jalview v 2.11 software (18).

MHC molecules are essential for presenting the pathogen-derived peptides to the respective T cells. Thus, it is necessary to perform molecular docking between T-cell epitopes and MHC molecules to understand the possible binding mode between them. First, the PEP-FOLD 3.5 server (https://bioserv.rpbs.univ-paris-diderot.fr/services/PEP-FOLD3/) (35)was used to model the three-dimensional (3D) structure of T-cell epitopes as the ligand, and RCSB Protein Data Bank (RCSB PDB) website (https://www2.rcsb.org/search/advanced/sequence) (36) was employed to acquire the 3D structure of MHC molecules as for the receptor. For some class I and class II MHC molecules for which the 3D structures are unavailable on the RCSB PDB website, we use the 3D structures of HLA-A*02:01 and HLA-DRB1*01:01 as the substitutes, respectively. Then, the energy minimization of receptors and ligands was carried out by Swiss-PdbViewer PDB viewer software (37), and following that, the molecular docking was run by Autodock Vina software hosted in PyRx v0.8 (38). AutoDock Vina, a classic molecular docking tool, has better performance for the prediction of docking between small molecules and proteins (39). The docking results were further analyzed using LigPlot+ v.2.2 software (40). Furthermore, to comparatively evaluate of binding affinity between the selected epitopes and MHC molecules, the molecular dockings between seven experimental verified antigen epitopes of poxvirus and HLA-A*02:01 molecules were selected as the positive controls.

All T-cell and B-cell epitopes that satisfied the filtration criteria were incorporated into the construction of the multi-epitope vaccine. The linker is an essential part of recombinant fusion proteins, which is capable of forming effective epitope conjugation between epitopes and allowing the independent immunological activities of the epitopes (41). Low immunogenicity is a major deficiency of multi-epitope vaccines. According to the previous study, a suitable adjuvant molecule can significantly boost the peptide vaccine’s immunogenicity and longevity (42). 50s ribosomal protein L7/L12 from Mycobacterium tuberculosis is capable of promoting DC maturation and diverse pro-inflammatory cytokine production via interacting with TLR4, which helps enhance cellular immunity (43). In addition to direct antiviral activity, β-defensins can regulate adaptive immunity to viral infection by recruiting naive T cells and immature dendritic cells (DCs) to the infected site (44). LT-IIC, a type II heat-intolerant enterotoxin (HLT) produced by Escherichia coli (E.coli), is an attractive adjuvant with the unique property of enhancing the body’s humoral and cellular immune responses to co-administered antigens (45). The Cholera Toxin B subunit (CTB) is also considered an attractive adjuvant due to its ability to efficiently bind to antigen-presenting cells and enhance the stability and half-life of the entire vaccine molecule (46). RS09, an LPS peptide mimic, is an agonist of TLR4. When used as an adjuvant in certain vaccines, it is capable of resulting in stronger immune activation and higher antibody production (47). All five adjuvants have been widely used in the design of multi-epitope vaccines in silico (48–51) and have been shown to stimulate the host’s immune response in animal experiments (43, 45, 52–54). So, in this study, each of the five adjuvants was added to the N-terminus of the vaccine using EAAAK liners, resulting in the formation of five vaccine constructs. Afterward, the PADRE sequence was linked to the adjuvant with the help of the EAAAK linker. When added to the vaccine construct, the PADRE, which can bind to most common HLA-DR molecules with high affinity, has been shown to optimize antibody responses and deliver significant helper T-cell activity in vivo ( 55). The CTL epitopes were connected by AAY linkers, the HTL epitopes were linked using GPGPG linkers, and the B-cell epitopes were spaced with the help of KK linkers. The last CTL epitope was linked to the first HTL epitope and the first HTL epitope was connected to the first B-cell-epitopes using the HEYGAEALERAG linker. Finally, the “6xHis tag” was appended to the C-terminal of the vaccine construct via the RVRR linker to facilitate protein purification but not affect the functionality of the fusion protein.

Out of all vaccine constructs, only antigenic, non-allergenic, non-toxic, stable, hydrophilic, thermostable, without transmembrane domains and signal peptide sequence, and highly soluble constructs were further submitted to the PSIPRED 4.0 server (http://bioinf.cs.ucl.ac.uk/psipred/) to obtain its secondary structure. The PSIPRED 4.0 server performed protein secondary structure prediction based on position-specific scoring matrices with a prediction accuracy of over 84% (60). The tertiary structures of the vaccine constructs were modeled using the Robetta server (http://robetta.bakerlab.org). Based on the fact that a confident template for the putative domains of a submitted sequence is available or unavailable, this server performs comparative modeling or de novo modeling, respectively (61). And this server will generate five models of each submitted sequence. Then, the optimal crude 3D models of each vaccine sequence were submitted to the GalaxyRefine module placed in the GalaxyWEB server(https://galaxy.seoklab.org/), which can rebuild unreliable loops or termini of the initial model structures using an optimization-based refinement method and further generate five refined models of each crude model (62). Referring to the criterion that the better model with a lower MolProbity score, we picked the best 3D refined model for each vaccine construct and submitted it to the ERRAT tool, PROCHECK tool of structure validation service SAVES v6.0 (https://saves.mbi.ucla.edu/) and ProSA-Web server (https://prosa.services.came.sbg.ac.at/prosa.php) for quality verification. The overall quality factor generated by the ERRAT program is the number of non-bonded interactions (side chains) formed between pairs of different atomic types within 0.35nm, and models with an overall quality factor value greater than 85 are generally considered to have good quality (63). The Ramachandran plot generated by the PROCHECK tool checks the stereochemical quality of the protein structure, and a better model has more residues in the Ramachandran-favored region and fewer residues in the Ramachandran-disallowed region. The Z-score predicted by the ProSA-Web server reflects the overall quality of the model, and a positive Z-score indicates the input structure has problematic or erroneous parts (64).

Stable interaction between the vaccine candidate and immune receptor is an important precondition to generating successful immune reactions. The TLRs not only hold a vital position in the first line of defense against pathogens but also play a vital role in linking innate immunity with adaptive immunity (65). Previous studies revealed that TLR2 and TLR4 can induce antiviral immune responses by detecting the virus coat proteins (66, 67). Class I and class II MHC molecules help the CD4+ T and CD8+ T cells recognize foreign antigens(vaccine). Therefore, the molecular docking of the vaccine candidates with TLRs (TLR2, TLR4) and MHCs (HLA-A*02:01, HLA-DRB1*01:01) was performed using the ClusPro 2.0 server (https://cluspro.bu.edu/login.php) (68). Since 2004, ClusPro has consistently been one of the best docking servers as demonstrated in Critical Assessment of Predicted Interactions (CAPRI), providing high predictive performance for the docking of protein-protein complexes (68). In each docking case, a total of thirty models were generated, and ten models with highly populated clusters were presented. Then, the best model in each docking case was further submitted to the HADDOCK 2.4 server (https://wenmr.science.uu.nl/haddock2.4/) (69), which is capable of refining the complex structures by rearranging the side-chain in the restricted interface and optimizing soft rigid-body, for refinement. After refinement, the docked models were submitted to the PRODIGY server (https://bianca.science.uu.nl/prodigy/) (70) for predicting binding energy. Furthermore, the PDBsum server (http://www.ebi.ac.uk/thornton-srv/databases/cgi-bin/pdbsum/GetPage.pl?pdbcode=index.html) (71) was utilized to investigate interacting residues between docked chains (i.e., the vaccines, the TLRs and MHCs). The PyMOL v2.5 was used to visualize the vaccine-receptor complex (72), and the Blender v3.3 software was employed to render the final image (73).

More than a thousand distinct HLA alleles have been identified in humans, and the frequency of distribution of each HLA allele varies across regions and ethnicities. Given that T lymphocytes can only recognize molecular complexes composed of epitopes bound to MHC molecules, a particular epitope triggers immune responses only in individuals expressing HLA alleles with an affinity for this epitope. The IEDB population coverage tool (http://tools.iedb.org/population/) was employed to calculate the population coverage of the designed vaccine in this work based on the data of TCR specificity, HLA restriction, and HLA allele frequencies.

For assessing the potential immune efficacy of the designed vaccine, we submitted it to the C-IMMSIM v10.1 tool (https://kraken.iac.rm.cnr.it/C-IMMSIM/) (82). The C-IMMSIM server recognizes epitopes by Position Specific Scoring Matrix (PSSM) methods, and can simultaneously stimulate the immune response of three compartments that represent three separate anatomical regions found in mammals (the bone marrow, the thymus and lymph nodes) to antigen epitopes. During the immune simulation, all parameters except for time steps were left at their default settings. In this study, the simulation time steps were set to 1050 (1-time step equals 8 h), and each injection point was set at time steps 1, 84, and 168, respectively. The default dose of each antigen injection was 1,000, and the time interval between each dose was four weeks, which is the recommended interval between injections for most commercial vaccines. To comparatively evaluate the candidate vaccine’s immunogenicity, a novel multi-epitope vaccine against Echinococcus granulosus, whose antigenicity and immunogenicity had been verified in vitro and in vivo experiments (83), was chosen as a positive control.

Due to the degeneracy of codons, each amino acid may correspond to multiple codons. However, there are great differences in codon preference among different species and organisms. To acquire the codon that is capable of efficiently encoding the targeted amino acid in the selected expression host, the codon optimization of the vaccine was conducted using the JCat server (http://www.jcat.de/CAICalculation.jsp) (84). The JCat server is an easy-to-use program, which can perform real-time calculations, eliminate the need for the definition of highly expressed genes, and avoid cleavage sites for specific restriction enzymes (84). The Escherichia coli (E. coli) K12 strain was selected as the expression host. Two parameters, codon adaptation index (CAI), and the percentage of GC content, were used to evaluate the expression efficiency of the vaccine in E. coli. The CAI, with a value range of 0 to 1, refers to the coincidence degree between synonymous codons and optimal codons, and a higher CAI indicates a higher expression level of the foreign gene in the host. The percentage of GC content represents the guanine (G)+cytosine (C) ratio in the DNA sequence, and its optimal range is 30%–70%. Otherwise, undesirable gene expression levels and transcriptional efficiency will result. Afterward, the BamHI and XhoI restriction endonuclease sites were integrated into the N- and C-terminal sites of optimized codon sequences, respectively. The entire sequence was cloned into the pET28a (+) vector using SnapGene software. Finally, the MPXV-5 vaccine protein was synthesized by Sangon Biotech co ltd (Shanghai, China). The flowchart of vaccine design is illustrated in Figure 1 .

The NCBI accession numbers and amino acid sequences of A35R, B6R and H3L derived from three objective MPXV strains were presented in Supplementary Table 1 . Protein sequence alignment results ( Figure 2 ) demonstrated that there was no difference in three target protein sequences (A35R, B6R and H3L) between MPXV-UK_P2 and 2022 reference strain. Compared with the Zaire-96-I-16 MPXV strain, there were 3 mutation sites in A35R and 2 mutation sites in H3L, but no mutation sites in B6R. As shown in Figure 2 , there were two mutation sites falling into the epitope regions.

The immune system is comprised of many specialized cell types, all of which work together to keep people healthy. CTLs, the main force of cellular immunity, are immune cells that can directly mediate the death of target cells and eliminate circulating antigens in the host body (85). Stimulated B cells can differentiate into antibody-producing plasma cells that play a vital role in the generation of humoral immune responses (86). HTLs, another group of immune cells, can assist both humoral and cellular immunity by secreting different cytokines (87). Hence, to construct an effective multi-epitope vaccine, both T-cell epitopes (CTL and HTL epitopes) and B-cell epitopes should be incorporated. After the layer-by-layer screening ( Figure 3 ), a total of 7 CTL epitopes ( Table 1 ), 6 HTL epitopes ( Table 2 ), and 7 liner B-cell epitopes ( Table 3 ) were finally included in the vaccine constructs. The binding affinity (IC50) of the experimentally determined poxvirus-specific T cell epitopes to HLA molecules was shown in Supplementary Table 2 .

Except for only one epitope (STETSFNDK), all other epitopes had multiple binding alleles. Some had even as high as 8 (TMSAFLIVR) or 13 (FTYTGGYDV) alleles ( Table 1 ). In this study, we performed molecular docking between each T-cell epitope and one of their corresponding alleles for a representative docking analysis ( Supplementary Table 3 ). As shown in Supplementary Table 3 , Autodock Vina docking data indicated that the binding energy between CTL epitopes and MHC class I alleles ranged from -10.2 to -5.9 Kcal/mol, and the binding energy range between HTL epitopes and MHC class II alleles was -7 to -5.8 Kcal/mol, the binding energy range between seven experimental verified antigen epitopes of poxvirus and HLA-A*02:01 molecules was -9 to -5.8 Kcal/mol, indicating the T-cell epitopes had a good affinity for MHC molecules. The predicted results of Ligplus v2.25 demonstrated that many hydrogen bonds were formed between T-cell epitopes and MHC molecules, and salt bridges were also detected in certain docking cases ( Supplementary Table 3 ), which once again proved that T-cell epitopes had a strong affinity to MHC molecules.

The secondary structures of the final selected vaccine constructs were predicted by the PSIPRED server and predicted results ( Supplementary Figure 3 ) showed that the MPXV-2 contained 58.22% α-helix, 31.22% coils and 10.56% β-strands. Compared to MPXV-2, the MPXV-5 secondary structure comprised a greater proportion of α-helix (61.10%), a smaller proportion of β-strands (6.79%), with the similar proportion of coils (32.11%) ( Supplementary Figure 3 ).

The 3D structure of the selected vaccine constructs was predicted by the Robetta server and then refined by the GalaxyWEB server. Based on the model quality of the refined models ( Supplementary Table 5 ), the refined model 5 ( Figure 4A ) for MPXV-2 and the refined model 4 ( Figure 4A ) for MPXV-5 were chosen for further quality validation. The corresponding Overall Quality Factors generated by the ERRAT program were 91.8465 (MPXV-2-5) and 89.2105 (MPXV-5-4) ( Supplementary Figure 4 ), and the corresponding Z-scores predicted by the ProSA-Web server were -6.94 (MPXV-2-5) and -6.36 (MPXV-5-4) ( Figure 4B ). The Ramachandran plot analysis of the refined model MPXV-2-5 revealed that 92.3% of residues located in the most favorable region, 6.4% in the allowed region, and only 1.3% in the disallowed region ( Figure 4C ). For the refined model MPXV-5-4, 92.4% of residues suited in the most favored region, 6.4% in the allowed region, and only 1.2% of residues in the disallowed regions ( Figure 4C ). The above analysis confirmed the excellent quality of the refined 3D models of the MPXV-2 and MPXV-5 vaccine constructs.

The ClusPro generated thirty clusters in each docking case, and the best one with the largest cluster size in each docking case was selected to present ( Table 5 ). Table 5 displayed the detailed docking results for all refined docking complexes. Compared with the MPXV-2-immune receptor complexes, all corresponding MPXV-5-immune receptor complexes had lower docking scores and binding energy, indicating that MPXV-5 had stronger affinity for all immune receptors (TLR2, TLR4, HLA-A*02:01 and HLA-DRB1*01:01). The predicted results of the PDBsum server ( Figure 5 ; Supplementary Figure 5 and Table 5 ) showed that 6 hydrogen bonds, 11 hydrogen bonds and 2 salt bridges, 12 hydrogen bonds and 1 salt bridges, 16 hydrogen bonds and 6 salt bridges were formed in MPXV-2-TLR2, MPXV-2-TLR4, MPXV-2-HLA-A*02:01 and MPXV-2-HLA-DRB1*02:01 complexes, respectively. There are 12 hydrogen bonds with 3 salt bridges, 18 hydrogen bonds with 3 salt bridges, 12 hydrogen bonds with 3 salt bridges, and 12 hydrogen bonds with 5 salt bridges in MPXV-5-TLR2, MPXV-5-TLR4 and MPXV-5-HLA-A*02:01 and MPXV-HLA-DRB1*02:01 complexes, respectively ( Figure 5 ; Supplementary Figure 5 and Table 5 ). Docking studies showed that there were significant interactions between the designed vaccines (MPXV-2 and MPXV-5) and the immune receptors.

Based on the HLA alleles binding to T-cell epitopes in the multi-epitope vaccine, the population coverage predictions were performed to determine the proportion of the population in different parts of the world that would respond to the designed vaccines. As shown in Figure 8 , the population coverage of the multi-epitope vaccine was high in countries suffering relatively severe MPXV attacks in this outbreak, such as Europe (Spain, France, Germany, the United Kingdom, the Netherlands, Portugal), North America (Mexico, Canada, the United States), and South America (Brazil, Peru). The analyses of population coverage were suggestive of the fact that the designed vaccine candidates (MPXV-2 and MPXV-5) have the potential to combat the MPXV infection globally.

The immune simulation results of MPXV-2, MPXV-5 and control group (C1) were shown in Figure 9 ; Supplementary Figure 9 . Similar immune response patterns were observed between the simulations of candidate vaccines (MPXV-2 and MPXV-5) and C1. The primary response in the three simulations were both characterized by elevated immunoglobulin IgM ( Figure 9A ; Supplementary Figure 9A ), B cells (isotype: lgM) ( Figure 9B ; Supplementary Figure 9B ) and plasma cells (isotype: lgM) ( Supplementary Figure 9C ). In secondary and tertiary reactions, these elevations were elevated to a greater degree, and the elevation of different immunoglobulins and the immunocomplexes (IgM+IgG, IgG1+IgG2, IgG1, memory B cells and plasma cells (isotype: lgG1, lgG2) were also observed ( Figures 9B, C ; Supplementary Figure 9B, C ). As shown in Figures 9A–C ; Supplementary Figure 9A–C the titers of antibody (IgM + IgG, IgG1 + IgG2) and the level of plasma B cell (isotype: IgG1) in MPXV-2 group were higher than those in C1, while the titer of antibody (IgM + IgG, IgG1 + IgG2) and the level of plasma B cell (isotype: IgG1) in MPXV-5 group were lower than those in C1. Additionally, the helper T-cells along with corresponding memory cells were observed for two vaccine candidates and C1 ( Figure 9D ; Supplementary Figure 9D ). The levels of active cytotoxic T cells and resting cytotoxic T cells in MPXV-2 and MPXV-5 groups changed significantly with time, while the resting cytotoxic T cell level in C1 changed little, and no inactive cytotoxic t cells were observed in the C1 group ( Figures 9E ; Supplementary Figure 9E ). High levels of IFN-γ and IL-2 were observed in vaccine candidate groups and C1 group in the cytokines’ simulation plot ( Figure 9F ; Supplementary Figure 9F ), while the concentration level of IL-2 was slightly higher for MPXV-2. During the immune simulation, the increased level of NK cells, macrophage (MA) and DC cells were also detected in three groups ( Figure 9G–I ; Supplementary Figure 9G–I ). Together, the simulation results suggest that the MPXV-2 and MPXV-5 both can induce stronger immunity in the host body, thereby effectively combating the MPXV. Moreover, we can know that the humoral immune response provoked by MPXV-2 was much stronger than that of MPXV-5.

For the optimized codon sequences of MPXV-2 and MPXV-5 (Supplementary Data), the CAI were 0.97 and 0.94 respectively, and the GC content was 48.79% and 48.18%, respectively, showing that both the MPXV-2 and MPXV-5 might well express in E. coli K12 strain and have the potential for mass production. The optimized codon sequences of MPXV-2 and MPXV-5 were inserted into the pET28a (+) vector ( Supplementary Figure 10 ). Electrophoretogram of enzyme digestion of pET28a (+)-MPXV-5 was shown in ( Supplementary Figure 11 ).