Design and evaluation of a multi-epitope DNA vaccine against HPV16

ABSTRACT

Cervical cancer, among the deadliest cancers affecting women globally, primarily arises from persistent infection with high-risk human papillomavirus (HPV). To effectively combat persistent infection and prevent the progression of precancerous lesions into malignancy, a therapeutic HPV vaccine is under development. This study utilized an immunoinformatics approach to predict epitopes of cytotoxic T lymphocytes (CTLs) and helper T lymphocytes (HTLs) using the E6 and E7 oncoproteins of the HPV16 strain as target antigens. Subsequently, through meticulous selection of T-cell epitopes and other necessary elements, a multi-epitope vaccine was constructed, exhibiting good immunogenic, physicochemical, and structural characteristics. Furthermore, in silico simulations showed that the vaccine not only interacted well with toll-like receptors (TLR2/TLR3/TLR4), but also induced a strong innate and adaptive immune response characterized by elevated Th1-type cytokines, such as interferon-gamma (IFN-γ) and interleukin-2 (IL2). Additionally, our study investigated the effects of different immunization intervals on immune responses, aiming to optimize a time-efficient immunization program. In animal model experiments, the vaccine exhibited robust immunogenic, therapeutic, and prophylactic effects. Administered thrice, it consistently induced the expansion of specific CD4 and CD8 T cells, resulting in substantial cytokines release and increased proliferation of memory T cell subsets in splenic cells. Overall, our findings support the potential of this multi-epitope vaccine in combating HPV16 infection and signify its candidacy for future HPV vaccine development.

Highlights

• Through the stringent selection of T-cell epitopes and other necessary elements, a novel multi-epitope vaccine targeting HPV 16 E6 and E7 oncoproteins was constructed using an immunoinformatics approach.

• The vaccine designed can induce both cellular and humoral immune responses, encompassing all the required immunogenic, physicochemical, and structural characteristics for an ideal vaccine design. Moreover, it offers decent worldwide coverage.

• In animal studies, the vaccine demonstrated strong immune responses, including expansion of CD4 and CD8 T cells, cytokine release, and enhanced memory T cell proliferation, resulting in long-term anti-tumor effects, inhibition of tumor growth, and prolonged survival in tumor-bearing mice.

• The immunological evaluation of the designed vaccine suggests its potential as a novel vaccine candidate against HPV 16.

GRAPHICAL ABSTRACT

KEYWORDS:

• HPV
• multi-epitope vaccine
• cancer vaccine
• immunoinformatics
• immune simulation

Introduction

Cervical cancer, ranked as the fourth most prevalent and fatal cancer in women globally, poses a significant public health challenge.¹ HPV is the primary culprit, with over 200 genotypes classified into low-risk (e.g., HPV6, 11, 40) and high-risk strains (e.g., HPV16, 18, 31, 33), the latter being implicated in various malignancies, particularly cervical cancer. HPV16, identified as the most potent oncogenic type, accounts for almost half of all cervical cancer cases.² Persistent infection with high-risk HPV strains is a crucial factor in the development of cervical cancer. While the prophylactic HPV vaccine effectively prevents viral infection, it does not eliminate existing infections.³ In contrast, therapeutic HPV vaccines aim to trigger T-cell-mediated immune responses in those already infected, potentially halting or reversing the disease’s progression.⁴

The continuous expression of HPV proteins E6 and E7 in infected cells is crucial for the development and maintenance of HPV-associated malignancies, particularly cervical cancer. These proteins play a significant role in the multistage transformation process of cervical cells by interacting with various cellular signaling elements that regulate the cell cycle, genomic stability, and epigenetic modifications.⁵ Specifically, E6 inhibits the protein p53, while E7 inhibits the retinoblastoma protein (pRb), leading to interference with cell cycle checkpoint regulation. Additionally, E6 and E7 proteins, being exogenous, avoid issues related to immunological tolerance.⁶ Therefore, the present therapeutic HPV vaccines target both the E6 and E7 proteins⁷ as well as either E6⁸ or E7⁹ individually.

Effective immune responses against tumors rely on the recognition and activation of CTLs in response to tumor antigens.^10,11 CTL epitopes, short peptides within antigen molecules, are recognized by the T-cell receptor on CD8+ T cells, leading to CTL production and immune response activation. These CTLs act as main effector cells, mediating the destruction of tumor cells by inducing cytotoxic responses and tumor cell apoptosis. HTL epitopes also play a crucial role in generating and maintaining antitumor immune responses by activating Th cells and enhancing CTL responses.^12–14 Immunoinformatics has facilitated the discovery and utilization of tumor antigenic T-cell epitopes for therapeutic tumor vaccine design. This approach, successfully applied in developing multi-epitope vaccines against various pathogens, including Hepatitis C virus,¹⁵ SARS-CoV-2,¹⁶ and Mycobacterium tuberculosis.¹⁷ Unlike conventional vaccine discovery methods involving labor-intensive pathogen cultivation, the emerging reverse vaccinology strategy, rooted in immunoinformatics, offers a promising, efficient approach for vaccine development.¹⁸

This study employed a computer network-based prediction method to identify T-cell immunogenic epitopes from HPV16 E6 and E7 oncoproteins. These epitopes were connected with suitable linkers and an adjuvant to create a multi-epitope chimera with potential resistance against cervical cancer. The physicochemical properties, antigenicity, allergenicity, toxicity, and solubility of the multi-epitope vaccine were thoroughly assessed. Immune simulation, computational cloning, and normal mode analysis evaluated the vaccine’s stability and effectiveness. Finally, the effectiveness of the vaccine was verified in animal experiments. Integrating reverse vaccinology, immunoinformatics, and computational simulations, this methodology was validated through mouse experiments, showcasing its potential for the development of an innovative multi-epitope chimeric vaccine targeting cervical cancer.

Materials and methods

Protein sequence retrieval

The amino acid sequences of the E6 protein (NP_041325.1), E7 protein (NP_041326.1), and β-defensin 3 (1KJ6_A) were obtained from the National Center for Biotechnology Information database (https://www.ncbi.nlm.nih.gov/).

Epitope prediction and population coverage analysis

The Immune Epitope Database (IEDB) was utilized to predict MHC-restricted T-cell epitopes from HPV16 E6 and E7 protein sequences. CTL epitopes for HLA class I alleles (HLA-A*02:01, HLA-A*24:02, HLA-A*11:01) were predicted using the “MHC-I Processing Predictions” and “MHC-I Binding Predictions” modules, employing the NetMHCpan and NetMHCpan EL 4.1 methods with specified thresholds. The NetMHCpan¹⁹ method was used for “MHC-I Processing Predictions” with thresholds set at −1.5, while the NetMHCpan EL 4.1²⁰ method was applied for “MHC-I Binding Predictions” with a threshold of 2. For HTL epitopes, the IEDB’s “MHC-II Binding Predictions” module was used to predict epitopes with high binding affinity to human HLA class II alleles (HLA-DRB1*03:01, HLA-DRB1*07:01, HLA-DRB1*15:01, HLA-DRB3*01:01, HLA-DRB3*02:02, HLA-DRB4*01:01, and HLA-DRB5*01:01).²¹ Similarly, predictions were made for C57BL/6 mouse MHC molecules (H2-Db, H2-Kb, H2-IAb), and antigenic peptides presented by both human and C57BL/6 mouse MHC molecules were identified.

The antigenicity of CTL epitopes was analyzed using the VaxiJen server²² (http://www.ddg-pharmfac.net/vaxijen/VaxiJen/VaxiJen.html); the potential of HTL epitopes to induce IFN-γ was assessed using the IFNEpitop platform²³ (http://crdd.osdd.net/raghava/ifnepitope/); and the toxicity of all epitopes was evaluated in the ToxinPred server²⁴ (http://crdd.osdd.net/raghava/toxinpred). The “Population Coverage” tool²⁵ (http://tools.iedb.org/population/) was employed to conduct population coverage analysis. CTL and HTL epitopes, utilized in vaccine construction, were individually and collectively analyzed in terms of population coverage.

Molecular interaction analysis of selected CTL and HTL epitopes with HLA alleles

The I-TASSER server²⁶ was utilized to predict the three-dimensional (3D) models of the selected epitopes. Furthermore, we obtained the crystal structures of HLA-A*02:01 (PDB ID: 1DUZ) and HLA-DRB1*01:01 (PDB ID: 1AQD) from the Protein Data Bank (PDB) (https://www.rcsb.org/). Subsequently, the selected CTL and HTL epitopes were docked to HLA-A*02:01 and HLA-DRB1*01:01, respectively, using the ClusPro 2.0 server²⁷ (https://cluspro.org/login.php).

Multi-epitope vaccine design and immunological and physiochemical properties

Then, we independently analyzed all candidate CTL epitopes to assess their compatibility for efficient binding to other epitopes to determine the order of epitopes in the final vaccine. A spacer sequence AAY²⁸ was added to the C-terminus of each CTL epitope, and I-TASSER predicted their 3D structures. The HADDOCK 2.4 server²⁹ (https://bianca.science.uu.nl/haddock2.4/) determined the compatibility between epitopes. The cluster with the maximum interaction compatibility of two epitopes was refined for compatibility with the third, and so on iteratively until the final vaccine structure was established. The spacer sequence GPGPG was introduced between the CTL and HTL epitopes and between the HTL epitopes.³⁰ The β-defensin 3 served as an adjuvant, linked ahead of the first CTL epitope via the EAAAK linker,³¹ forming the vaccine sequence named BD3-12P.

Physicochemical properties, antigenicity, allergenicity, toxicity, and solubility of HSP70-12P were predicted using the ProtParam server³² (https://web.expasy.org/protparam/), VaxiJen server,²² AllerTOPv.2.0 server³³ (https://www.ddg-pharmfac.net/AllerTOP/method.html), ToxinPred server,²⁴ and Protein-Sol server³⁴ (https://protein-sol.manchester.ac.uk/), respectively.

Structures modeling, validation, and refinement

PDBsum³⁵ (http://www.ebi.ac.uk/thornton-srv/databases/cgi-bin/pdbsum/GetPage.pl?pdbcode=index.html) and AlphaFold2 tool³⁶ were utilized to predict the Secondary structure and 3D models of BD3-12P, respectively. Subsequently, the resulting 3D model was refined using the GalaxyRefine server.³⁷ Structural validation was conducted using the PDBsum graphical database³⁸ (http://www.ebi.ac.uk/thornton-srv/databases/pdbsum/Generate.html) for Ramachandran diagrams, and the overall quality score was assessed via the ProSA-web server.³⁹

Molecular docking and simulation

The ClusPro 2.0 server²⁷ was utilized for analysis of the interaction of BD3-12P with the immune receptors TLR2 (PDB ID: 2Z7X), TLR3 (PDB ID: 2A0Z), and TLR4 (PDB ID: 4G8A). The PyMol visualization program was used to analyze the docking results. We conducted an analysis of the docked complexes involving the BD3-12P vaccine and TLR2, TLR3, and TLR4 receptors using the iMODS tool⁴⁰ (https://imods.iqfr.csic.es/). The C-ImmSim simulation server⁴¹ (http://150.146.2.1/C-IMMSIM/index.php) was utilized to simulate changes in the expression of T and B lymphocytes, as well as cytokines, after stimulation with BD3-12P. The workflow followed in this study is depicted in Supplementary Figure 1.

DNA vaccine preparation and experimental models

In order to make the vaccine sequence better expressed in eukaryotes, we sequentially added the Kozak sequence containing the start codon (GCCACCATG) and a piece of IgE leader sequence (DWTWILFLVAAATRVHS) at its 5’ end. The JCat codon optimization tool⁴² (http://www.jcat.de/) was utilized for reverse translation and codon optimization of BD3-12P. The codon-optimized sequence synthesized by Shanghai Generay Biotech, was inserted into the pVAX1 plasmid (Invitrogen, 35–1049) between the EcoRI and XbaI sites.

Female C57BL/6 mice were procured and housed in a pathogen-free environment at the Institute of Medical Biology, Chinese Academy of Medical Sciences. The Institute of Medical Biology’s Institutional Animal Ethics Committee approved all animal-based experiments (Approval number: DWSP202306014). TC-1 tumor cells, originating from C57BL/6 mice primary lung epithelial cells and transformed with c-Ha-ras and HPV-16 E6- and E7-encoding genes,⁴³ served as a preclinical tumor model for assessing therapeutic HPV vaccines targeting HPV16 E6/E7 antigens.

Immunization and tumor challenge

Mice were vaccinated via intradermal microneedle injection route. In all cases, mice were given three doses sequential vaccinations 7 days apart, each injection is 80 μg. In therapeutic (n = 7 per group) and survival experiments (n = 10 per group), mice were subcutaneously injected into the left flanks with TC-1 cells (1 × 10⁵ cells per mouse). When the TC-1 Challenge Day 7, different treatments were given to each group. Tumor growth was monitored, and volume was measured every 3–5 days with calipers. The tumor volume was calculated with the formula (L x W²/2), where L is the length and W is the width of the tumor. Mice were euthanized when the tumor volume reached 1500 mm³.

Flow cytometry

Spleen single-cell suspensions underwent flow cytometry staining for cell analysis. The following monoclonal antibodies were used for extracellular staining: CD3 (Invitrogen, 363-0032-82), CD4 (BioLegend 100,406), CD8 (BioLegend 100,752), NK1.1 (BioLegend 108,761), CD62L (BioLegend 104,412), CD44 (BioLegend 103,007). Intracellular staining utilized antibodies for IL2 (BioLegend 503,808), IL17A (BioLegend 506,926), IL4 (BD Biosciences 554,436), granzymeB (BioLegend 372,204), IFN-γ (BioLegend 505,808), and TNF-α (BioLegend 506,328). Stimulation and Cytokine Staining: Spleen cell suspensions (5 × 10⁶ cells per well) were stimulated with a peptide pool composed of 12 short peptides (10 µg mL⁻¹/peptide) for 4 hours. Subsequently, Brefeldin A (BioLegend 420,601) was added for overnight incubation. For surface marker staining, nonspecific signals were blocked using TruStain FcXM (BioLegend 101,320) on ice for 15 minutes. Viability staining followed the instructions of the Zombie NIR™ Fixable Viability Kit (BioLegend 423,106). Subsequently, surface markers were stained. For intracellular cytokine staining, fixation and permeabilization were performed, followed by incubation with cytokine-specific antibodies. Flow cytometry data were acquired on a Beckman Coulter flow cytometer, and analyzed using FlowJo 10.8.1 software (https://www.flowjo.com/solutions/flowjo/downloads).

Elispot assay

The Mouse IFN-γ pre-coated ELISpot plate (Product code: 3321-4AST-10) was washed, incubated with RPMI-1640 medium, and then loaded with spleen cells (2 × 10⁵ cells per well) and 12 short peptides (5 µg/mL per peptide). After 40 hours of incubation at 37°C in a 5% CO₂ environment, spot formation units of IFN-γ secretion were detected using the corresponding antibodies and substrate solution. Spot counts were obtained with an automated ELISpot image analyzer (CTL, S6 FluoroCore), and results were expressed as spot-forming units per 1 × 10⁶ cells after subtracting negative control values.

Statistical analysis

Statistical analyses utilized GraphPad Prism v.9.0 software (https://www.graphpad.com/). Unpaired two-tailed t-tests were employed to compare control and vaccine groups, with survival rates assessed via log-rank (Mantel-Cox) tests depicted in survival curves. Values and error bars denote means ± SD. Significance levels were defined as p < .05 (*p < .05, **p < .01, ***p < .001, and ****p < .0001).

Results

Identification and analysis of T-cell epitopes

A total of 10 CTL epitopes and 2 HTL epitopes were identified, with 98.36% and 78.51% coverage of the global population, respectively (Supplementary Table 1). Moreover, all epitopes were predicted to be nontoxic, 70% of CTL epitopes were identified as “Probable ANTIGEN” (with a threshold of 0.4, scores greater than 0.4 were considered as “Probable ANTIGEN”), and all HTL epitopes were identified as having the potential for IFN-γ induction (Table 1). The AHYNIVTF epitope showed the most favorable docking energy score of −953.5 (kcal/mol) with HLA-A*02:01 among the CTL epitopes. Similarly, among the HTL epitopes, the DLCIVYRDGNPYAVC and HEYMLDLQPETTDLY epitopes exhibited the lowest docking energy scores of −813.2 (kcal/mol) and −868.0 (kcal/mol), respectively, when docked with the HLA-DRB1*01:01 molecule (Supplementary Table 2).

Table 1. Overall the HTL and CTL epitopes selected to construct the vaccine candidate.

Protein	Name	Start	End	Length	CTL epitope	Alleles	Population coverage (World)	Reference	Assay description in the References	Antigenicity	Toxicity
E6	E61	86	95	10	YCYSLYGTTL	HLA-A*24:02, HLA-A*23:01	26.18%	/	/	0.4484	Non-Toxin
	E62	81	90	10	SEYRHYCYSL	HLA-B*40:01, HLA-B*44:02, HLA-B*44:03, HLA-A*23:01, HLA-A*24:02, HLA-B*08:01, HLA-A*30:01	50.83%	【46】	ICS TNF-α release Positive, ICS IFN-γ release Positive	1.7034	Non-Toxin
	E63	75	83	9	KFYSKISEY	HLA-A*30:02, HLA-A*30:01, HLA-A*03:01, HLA-B*15:01, HLA-A*32:01, HLA-A*31:01, HLA-A*23:01, HLA-A*24:02, HLA-A*11:01, HLA-B*35:01, HLA-B*58:01, HLA-A*68:01, HLA-A*33:01, HLA-A*26:01, HLA-B*57:01, HLA-B*44:03, HLA-B*44:02, HLA-A*01:01, HLA-B*08:01, HLA-A*02:03, HLA-A*02:06	89.32%	【47】	ICS IFN-γ release Positive	0.2257	Non-Toxin
	E64	124	132	9	RHLDKKQRF	HLA-A*23:01, HLA-A*24:02, HLA-A*30:02	28.27%	/	/	0.3547	Non-Toxin
	E65	49	57	9	VYDFAFRDL	HLA-A*23:01, HLA-A*24:02	26.18%	【48, 49】	51 chromium cytotoxicity Positive, ELISA IFN-γ release Positive, ELISPOT IFN-γ release Positive, ICS IFN-γ release Positive, multimer/tetramer qualitative binding Positive	1.4227	Non-Toxin
	E66	38	45	8	VYCKQQLL	HLA-A*23:01, HLA-A*24:02	26.18%	/	/	0.5081	Non-Toxin
E7	E71	7	15	9	TLHEYMLDL	HLA-A*02:03, HLA-A*02:01, HLA-A*02:06	41.35%	【50, 51】	ICS IFN-γ release Positive, in vitro assay cytotoxicity Positive, in vivo assay decreased disease Positive-Low, multimer/tetramer qualitative binding Positive	0.3274	Non-Toxin
	E72	77	87	11	RTLEDLLMGTL	HLA-B*40:01, HLA-B*58:01, HLA-B*57:01, HLA-A*02:06, HLA-A*32:01, HLA-A*02:01, HLA-A*30:01, HLA-A*02:03	55.46%	【52】	ELISPOT IFN-γ release Positive	0.5722	Non-Toxin
	E73	49	57	9	RAHYNIVTF	HLA-A*32:01, HLA-B*15:01, HLA-B*58:01, HLA-B*57:01, HLA-B*35:01, HLA-A*30:01, HLA-A*23:01, HLA-A*24:02, HLA-A*30:02, HLA-B*53:01, HLA-B*07:02, HLA-B*08:01, HLA-A*02:06, HLA-B*40:01, HLA-A*31:01, HLA-A*03:01, HLA-B*44:03, HLA-B*51:01	81.45%	【50, 53–55】	biological activity degranulation Positive, ELISA IFN-γ release Positive, ELISPOT IFN-γ release Positive, ICS IFN-γ release Positive, in vitro assay cytotoxicity Positive, in vitro assay proliferation Positive, in vivo assay decreased disease Positive, in vivo assay survival from pathogen challenge Positive, in vivo assay tumor burden after challenge Positive, multimer/tetramer qualitative binding Positive	0.5919	Non-Toxin
	E74	50	57	8	AHYNIVTF	HLA-A*23:01, HLA-A*24:02, HLA-B*15:01, HLA-B*44:03, HLA-B*35:01, HLA-A*32:01, HLA-A*30:02, HLA-B*44:02, HLA-B*58:01	53.73%	/	/	0.5177	Non-Toxin
Protein	Name	Start	End	Length	HTL epitope	Alleles	Population coverage (World)	Reference	Assay description in the References	IFN-γ inducing	Toxicity
E6	E67	56	70	15	DLCIVYRDGNPYAVC	HLA-DRB3*01:01, HLA-DRB1*03:01, HLA-DRB3*02:02, HLA-DRB1*13:02, HLA-DRB1*04:01	33.64%	/	/	POSITIVE	Non-Toxin
E7	E75	9	23	15	HEYMLDLQPETTDLY	HLA-DRB3*01:01, HLA-DQA1*05:01/DQB1*02:01, HLA-DQA1*03:01/DQB1*03:02, HLA-DQA1*01:01/DQB1*05:01	67.61%	/	/	POSITIVE	Non-Toxin

Multi-epitope vaccine design and immunological and physiochemical properties

The final sequence of CTL epitope connections obtained through this comprehensive combination analysis was E61-E64-E74-E71-E73-E62-E72-E66-E65-E63 (Supplementary Table 3). The selected epitopes and adjuvant were concatenated in the order obtained from the above analysis using EAAAK, AAY, and GPGPG linkers (Figure 1a). BD3-12P had a practical molecular weight and physiochemical properties. The aliphatic index was high (75.79), possibly indicating the protein’s thermal stability. The GRAVY score (−0.248) indicated the protein’s overall hydrophilic nature with better potential for interaction with neighboring water molecules. The solubility of BD3-12P was possibly suitable (0.49), surpassing the threshold value of 0.45. The final sequence demonstrated high antigenicity modeled on the virus (0.5620), and BD3-12P was predicted to be nontoxic and non-allergenic. (Supplementary Table 4).

Figure 1. Structural information diagram of BD3-12P. (a) Schematic representation of BD3-12P. (b) Predicted 2D structure. (c) The comparison between unrefined in cyan and refined in orange tertiary structures, and overlap of the two structures. (d) Evaluation of the 3D model before and after refining. The Z score for the unrefined model was calculated as − 1.31, which improved to − 1.79 in the refined model. The Ramachandran plot analysis shows that in the unrefined model, 77.0% of residues were in the favorable region, 17.6% in the additional permissive region, 2.7% in the generous permissive region, and 2.7% in the nonpermissive region. In the refined model, these values changed to 87.7%, 10.2%, 1.1%, and 1.1%, respectively.

Structures modeling, validation, and refinement

The results revealed that BD3-12P is composed of 6 helices, 3 helix-helix interactions, 25 β-turns, and 10 γ-turns (Figure 1b). The AlphaFold2 tool predicted five models for BD3-12P, and Model 1 was selected for structural improvement to have a tertiary structure much closer to the natural state. To assess the quality of BD3-12P and determine the differences between the unrefined and refined structures, we superimposed the two structures (Figure 1c). The refined model showed improvement compared to the unrefined model in terms of the Ramachandran plot analysis and overall quality score calculated using the ProSA-web server (Figure 1d).

Molecular docking and molecular dynamic simulation

The server generated 30 clusters per complex, ranked by maximum cluster size. TLR2, TLR3, and TLR4 interactions with BD3-12P were assessed, with cluster 0.00 having maximum sizes of 86, 53, and 76. Corresponding minimum energy scores were −1476.7, −1077.9, and −1306.0 Kcal/mol, respectively (Supplementary Table 5). PyMOL visualized docking results, and Ligplot + v.2.2.5 analyzed 2D interactions. The results indicate that BD3-12P interacts differently with TLR2, TLR3, and TLR4. Specifically, the A chain of TLR2 forms 4 hydrogen bonds and 2 salt bridges with BD3-12P, whereas the B chain establishes 14 hydrogen bonds and 1 salt bridge (Supplementary Figure 2). Additionally, the interaction of BD3-12P with the A chain of TLR3 results in the formation of 16 hydrogen bonds and 1 salt bridge (Supplementary Figure 3). Similarly, interactions with the A chain of TLR4 result in 19 hydrogen bonds, while no salt bridges are observed. Conversely, the B chain of TLR4 forms 11 hydrogen bonds and 3 stabilizing salt bridges (Figure 2).

Figure 2. Docked complex of BD3-12P and TLR4. (a) A Chain, B Chain, C Chain, and D Chain of TLR4 are shown in cyan, blue, magenta, and green, respectively, while BD3-12P is shown in red. (b-c) The map of interactions between BD3-12P and A Chain and B Chain of TLR4, in which, hydrogen bonds and salt bridge interactions were represented by dashed lines in spring green and red, respectively, while hydrophobic interactions are indicated by curved semicircles.

The deformability graph highlights regions with high deformability in the docking complexes, indicated by prominent peaks. Peaks were less pronounced across all complexes, suggesting low deformability of individual amino acid residues, indicating stability. Eigenvalues, correlating with energy needed for deformation, were 7.582751e-06, 1.611329e-05, and 1.256640e-05 for TLR2-BD3-12P, TLR3-BD3-12P, and TLR4-BD3-12P, respectively. Higher values imply more energy for deformation, indicating greater stability for TLR3-BD3-12P and TLR4-BD3-12P. The covariance map reveals correlated motions (red) in all complexes, indicating coordinated residue movements. B-factor values, reflecting internal vibrations, align with root-mean-square fluctuations, supporting complex stability. The elastic network model identifies atom pairs linked by springs, with darker regions indicating greater elasticity. Rigid regions in all complexes contribute to overall structure maintenance. In summary, these results suggested high stability in BD3-12P vaccine-TLR complexes during molecular dynamics simulations, with TLR3-BD3-12P and TLR4-BD3-12P showing stable internal vibrations and synergistic motions, enhancing our understanding of their dynamic behavior and potential biological roles (Figure 3).

Figure 3. The analysis of internal vibrational modes, conformational changes, and correlations within the docked complexes TLR2-BD3-12P (a), TLR3-BD3-12P (b), TLR4-BD3-12P (c).

Codon optimization and immune simulation

The optimized BD3-12P nucleotide sequence was cloned into the pVAX1 vector using SnapGene software (Supplementary Figure 4). Immune Simulations on C-ImmSim server explored four schedules (1, 2, 3, and 4 weeks). Results indicate that BD3-12P activates innate immune responses by stimulating macrophages and DCs (Supplementary Figure 5), and rapidly increases the numbers of active B cells, TH cells, and TC cells (Supplementary Figure 6).

BD3-12P induced strong humoral and cytokine immune responses, with antibody levels peaking after the third stimulation. Antigen levels peaked within two days post-injection, gradually declining as antibody production by B- and T-cells increased. The cytokine levels of IFN-γ, IL2, TGF-β, IL10, and IL12 also exhibited distinct peaks following each stimulation. IFN-γ levels rapidly rose after the initial immunization, reaching a consistent maximum value of 2.4 × 10⁶ ng/ml with 1- or 2-week intervals. Immunization at 3- or 4-week intervals resulted in a peak of 2.1 × 10⁶ ng/ml after the third stimulation, followed by a decline to near-zero levels one month later. The Simpson’s index consistently revealed a low D value after each injection, indicative of a robust diversity of BD3-12P (Figure 4). Based on these results, immunization intervals of 1 or 2 weeks may yield superior immune responses.

Figure 4. The levels of antibodies (IgG, IgM, IgG1, and IgG2) (a), cytokines and interleukins (b) after antigen 3 dose injection (black vertical lines) predicted by the C-ImmSim server.

Features of pVAX1-BD3-12P‐induced T‐cell immunity

To assess pVAX1-BD3-12P immunogenicity, three doses of the vaccine were administered to mice (Figure 5a). ELISpot analysis (Figure 5b) revealed robust IFN-γ secretion in splenocytes, indicating a potent response compared to the control. The functional role of CD4+ cells in the anti-tumor response was explored, showing substantial IL2, IFN-γ, and IL17A release from CD4 T cells in vaccinated mouse splenocytes after 12 HPV peptide mixture stimulation, while IL4 showed no significant change (Figure 5c). The impact on splenic CD8 T cell responses post-12-peptide mix stimulation demonstrated a higher frequency of CD8 T cells producing granzyme B, IFN-γ, TNF-α, or IL2 after triple vaccine immunization, indicating significant activation compared to the control group (Figure 5d).

Figure 5. Features of pVAX1-BD3-12P‐induced T‐cell immunity. (a) Experimental design. Mice (n = 5 per group) received 80 μg of pVAX1-BD3-12P intradermally with three weekly injections. (b) Evaluation of the T-cell immune response induced by the vaccine. The statistical analysis measured IFN-γ producing cell numbers, represented as spot-forming units (SFU). (c) Measurement of IL2, IFN-γ, IL4, and IL17A levels released by CD4+ T cells from splenocytes post-stimulation with the 12 HPV peptide mixture. (d) Measurement of IL2, IFN-γ, TNF-α, and granzyme B levels released by CD8+ T cells from splenocytes post-stimulation with the 12 HPV peptide mixture, assessed via flow cytometry. Data presented as mean ± standard deviation, analyzed using unpaired two-tailed t-tests (**P < .01, ***P < .001, ****P < .0001).

pVAX1-BD3-12P exhibits therapeutic efficacy and induced cellular immune response in a TC-1 murine model

To evaluate therapeutic efficacy, mice were challenged 1 × 10⁵ TC-1 cells on day 0 and intradermal immunizations with a 80 μg vaccine dose on day 7 (Figure 6a). ELISpot analysis demonstrated significantly induced IFN-γ secretion in pVAX1-BD3-12P-vaccinated mice compared to controls (Figure 6b). Tumor growth was markedly inhibited with pVAX1-BD3-12P, maintaining volumes below 100 mm³ by day 25, while control group tumors continued growing (Figure 6d). After the third immunization, mice were euthanized, and the spleen and tumor tissues were isolated. The results indicated that approximately 43% (3/7) of vaccinated mice showed tumor eradication, and others exhibited delayed growth, contrasting with larger tumors in the control group (Figure 6c).

Figure 6. pVAX1-BD3-12P mitigates or even eradicates tumor growth in the TC-1 mouse model. (a) Experimental design. Mice (n = 7 per group) received subcutaneous injections of TC-1 tumor cells (1 × 10⁵ cells per mouse). Mice harboring TC-1 tumors were immunized intradermally with 80 μg of pVAX1-BD3-12P on days 7, 14, and 21. (b) Evaluation of vaccine-induced T-cell immune response through ELISPOT, SFU denotes spot-forming units. Data are presented as mean ± standard deviation, and statistical significance was determined by unpaired two-tailed t-tests (**P< .01). (c) Representative images of tumor groups. (d) Kinetics of tumor growth. The green arrows indicate the time point at which the vaccine was injected. (e) Survival curve analysis of mice challenged with TC-1 tumor cells using the same immunization regimen and dosage as in figure a (n = 10 per group). Survival was analyzed using log-rank (Mantel-Cox) tests (****P < .0001).

Using the same immunization schedule and dose, a separate batch of mice showed extended survival in approximately 50% (5/10) of the experimental group beyond 60 days post tumor challenge, whereas all unvaccinated control group mice were euthanized within 35 days post TC-1 challenge due to extensive tumor growth (0/10) (Figure 6e).

Splenocytes, collected post-vaccination, underwent flow cytometry analysis to assess immune responses in tumor-bearing mice. Vaccinated mice exhibited a significant rise in CD4 and CD8 T cell numbers upon peptide pool stimulation compared to controls (Figure 7a,c). This increase correlated with heightened IL2, TNF-α, IFN-γ, and granzyme B production by both CD4 and CD8 T cells, indicating enhanced Th1 response and cytotoxicity. IL4 levels, however, only increased in CD4 T cells, not in CD8 T cells (Figure 7b,d). Additionally, vaccinated mice displayed a higher ratio of effector memory T cells in both CD4 and CD8 T cell subsets, suggesting differentiation from naive T cells due to repeated vaccination. No notable changes were observed in central memory T cells (Figure 7e). Moreover, vaccinated mice exhibited an increased proportion of NK cells compared to the control group (Figure 7f). Overall, pVAX1-BD3-12P vaccination in tumor-bearing mice elicited robust anti-tumor immune responses, promoting the proliferation of effector memory T cells within splenocytes. Figure 8 demonstrates the Gating Strategy for flow cytometry analysis.

Figure 7. pVAX1-BD3-12P stimulated the proliferation of CD4 T cells and CD8 T cells, along with their cytokine secretion in the spleens of TC-1 tumor-bearing mice, promoting the expansion of effector T cell subsets. (a) Evaluation of CD4 T cell proliferation in splenocytes. (b) Determination of IL2, IL4, TNF-α, IFN-γ, and granzyme B secretion in CD4 T cells from splenocytes. (c) Assessment of CD8 T cell proliferation in splenocytes. (d) Measurement of IL2, IL4, TNF-α, IFN-γ, and granzyme B secretion in CD8 T cells from splenocytes. (e) Collective analysis of naïve T cells, effector memory T cells, and central memory T cells in CD4 T cells and CD8 T cells. (f) Evaluation of NK cell proliferation in splenocytes. Data are presented as mean ± standard deviation, and statistical significance was determined using unpaired two-tailed t-tests, denoting significant differences as **P< .01, ***P< .001, and ****P< .0001.

Figure 8. Gating strategy for flow cytometry analysis. (a) Gating strategy for analyzing T cells in intracellular cytokine staining panel. (b) Gating strategy for populations of naïve T cells, effector memory T cells, and central memory T cells in CD4 T cells and CD8 T cells. Percentages indicate the percent of the parent population.

pVAX1-BD3-12P vaccination provides prophylactic efficacypVAX1

We initiated a protective experiment in C57BL/6 mice (Supplementary Figure 7a). Continuous monitoring of tumor size revealed complete absence of tumor growth in mice immunized with the vaccine, contrasting with rapid tumor growth in those immunized with pVAX1 (Supplementary Figure 7b). After 28 days of observation, significant tumor volumes in the control group mice necessitated euthanasia, underscoring that sustained HPV-specific immune responses persisted post-vaccination, effectively preventing TC-1 tumor development.

Discussion

The immunoinformatics approach, a crucial element in the reverse vaccinology strategy, revolutionizes vaccine discovery. Its efficiency lies in the rapid and accurate prediction and design of novel immunogens and vaccines. We potentially designed a candidate therapeutic vaccine against HPV16, based on predictions from immunoinformatics tools.

We adopted a comprehensive strategy, considering both ‘binding’ and ‘processing’ predictions. The ‘processing’ prediction encompassed proteasome processing, TAP transport, and MHC binding, providing an overall intrinsic potential score for each peptide as a T-cell epitope. We utilized NetMHCpan¹⁹ to assess MHC binding efficiency, considering binding stability in terms of ‘Binding Affinity,’ applicable to all MHC alleles and demonstrating excellent performance in recent comparisons.⁴⁴ For ‘binding’ prediction to MHC class I molecules, we employed NetMHCpan EL 4.1,²⁰ crucial for evaluating peptide-MHC complex stability by assessing the ability to ‘Elute Ligand.’⁴⁵ Population coverage was a critical factor, and predictions were based on alleles with high coverage, with the combined epitopes showing over 90% coverage in most regions. Notably, among the selected epitopes, six CTL epitopes (E62: SEYRHYCYSL,⁴⁶ E63: KFYSKISEY,⁴⁷ E65: VYDFAFRDL,^48,49 E71: TLHEYMLDL,^50,51 E72: RTLEDLLMGTL,⁵² E73: RAHYNIVTF^{50, 53–55} were previously reported, while four CTL epitopes and two HTL epitopes were newly identified in this study. In addition, the chosen CTL and HTL epitopes are compatible with both human and mouse MHC molecules. This unique prediction and selection approach allow for direct efficacy validation in tumor treatment assays using the C57BL/6 mice, eliminating the need for labor-intensive tasks like growing transgenic mice.

We obtained 12 T-cell epitopes through a unique prediction method, followed by further exploration of various strategies, including linking selected epitopes and designing candidate vaccines, to enhance the efficacy of our DNA vaccine, given reports in some studies suggesting potential low immunogenicity in DNA vaccines.^56–59 The BD3-12P formulated with multiple epitopes may demonstrate a strong tumor cell killing capability when expressed. Single-epitope targeting, while widely employed, is often limited in its ability to establish durable therapeutic effects. Conversely, multi-epitope vaccines, which are recombinant vaccines that address the limitations of single-epitope targeting, offer the advantage of generating more sustained epitope presentation and stronger immune responses. Furthermore, we incorporated HTL epitopes to synergistically enhance CTL immune responses by inducing activation of HTL cells, thus exerting a more robust and broader immune-killing effect on tumor cells.

In our study, we carefully arranged CTL epitopes, adjuvants, and linkers within the vaccine construct, ensuring optimal linkage. The final sequence of CTL epitope connections was obtained through a comprehensive compatibility analysis. In our approach, we utilized AAY linkers to fuse the CTL epitopes, which are known to be preferentially cleaved by the proteasome.²⁸ We utilized GPGPG linkers to fuse the HTL epitopes. This strategy aligns with the work of Living et al.³⁰ who designed a generic spacer sequence GPGPG based on the recognition of peptide motifs by the HLA-DR molecule and used it to link multiple HTL epitopes, minimizing junctional epitopes and facilitating efficient separation between epitopes. This connection is similar to the work of Waqas et al.⁶⁰ and Tarrahimofrad et al.⁶¹ In addition, the development of built-in adjuvants has become a key step in designing epitope vaccines to overcome obstacles such as enzymatic degradation and facilitate better recognition by the immune system of target cells.⁶² Fusion with potent adjuvants, such as Mycobacterium tuberculosis HSP, TLR ligands, herpes simplex virus glycoprotein D, β-defensins, has shown promising results in enhancing antigen-specific CD8 T cell responses.^63,64 In our study, we added β-defensin 3 to the N-terminal end of the final sequence using the EAAAK linker aims to activate human antigen-presenting cells.⁶⁵ The use of the EAAAK linker not only improved stability but also minimized attachment to other protein regions.³¹ Codon optimization enhances fusion protein expression to improve the immunogenicity of DNA vaccines, increasing survival rates in tumor-bearing mice.⁶⁶ In our study, codon optimization was employed. In addition, IgE leader sequences were added to enhance epitope presentation and expression. In addition to these strategies, effective epitope selection and careful vaccine design are crucial considerations.

After optimization with the aforementioned strategies, BD3-12P was predicted to be potentially stable, antigenic, possibly non-allergenic, and possibly nontoxic, making it a promising vaccine candidate. It has a suitable molecular weight and instability index. Its half-life in mammalian reticulocyte was determined to be 30 hours, longer than the construct designed in the study by Sarkar et al.⁶⁷ BD3-12P is also hydrophilic and capable of interacting with water molecules because the GRAVY value was −0.248. However, in the study by Negahdaripour et al.⁶⁸ with a GRAVY value of 0.252, indicating a need for micelles to increase vaccine interactions in a polar environment. Refinement of the 3D structural model significantly enhanced its quality, with 87.7% of residues in the favorable region. According to the results of Samira et al.⁶⁹ in the refined model, this value changed to 60.8%. Additionally, the Z score for the refined model (−1.79) further attests to the improved quality. The improved quality suggests increased efficacy and safety potential for the vaccine candidate.

The expression of the TLR4 immune receptor in human cervical cancer HeLa cells was found to be upregulated by 100 times compared to other TLRs, suggesting a correlation between TLR4 and cervical cancer progression, as well as resistance to apoptosis.⁷⁰ Meanwhile, studies have also indicated that TLR2 and TLR3 may recognize certain endogenous ligands, hinting at a potential role for TLR2 and TLR3 in tumor immune surveillance.⁷¹ Molecular docking analyses suggested that BD3-12P potentially exhibited stronger binding affinity for TLR2 and TLR4 compared to TLR3, supported by lower energy scores and larger cluster sizes. Two-dimensional interaction analyses demonstrated that the vaccine formed more hydrogen bonds and salt bridges with TLR2 and TLR4, indicating potentially stronger binding affinity. Salt bridges, in particular, play a pivotal role in enhancing the stability of protein-protein complexes. Normal mode analysis results underscored the high stability of the BD3-12P vaccine-TLR complex throughout the simulation, with a special emphasis on the TLR3-BD3-12P and TLR4-BD3-12P complexes, which exhibited more stable internal vibrational modes and coordinated motion. Higher eigenvalues signify a greater energy requirement for substantial deformation. Specifically, the eigenvalues for the TLR2-BD3-12P, TLR3-BD3-12P, and TLR4-BD3-12P complexes were measured as 7.582751e-06, 1.611329e-05, and 1.256640e-05, respectively. For comparison, the eigenvalues for the TLR2-C543P and TLR4-C543P complexes in the study by Cheng et al.¹⁷ were 8.659505e-06 and 9.469126e-06, respectively. The EBV vaccine-TLR4 complex in Victor Omoboyede et al.‘s study⁷² yielded an eigenvalue of 4.987694e-06. Additionally, PV1A and PV3B complexes with TLR2 and TLR4, as investigated by Manisha Pritam et al.⁷³ exhibited eigenvalues of 1.064e-06 and 7.498e-09, respectively. These results suggest that the TLR3-BD3-12P and TLR4-BD3-12P complexes possess notably stable properties.

Additionally, compared to previous studies, the levels of IFN-γ and IL2 induced by BD3-12P, as revealed by our immune simulation results, highlight its strong potential for immune stimulation. Depending on the different immune protocols, the maximum values of IFN-γ and IL2 were predicted to range from 2.1–2.4 × 10⁶ ng/ml and 2.2–2.9 × 10⁶ ng/ml, respectively. In contrast, in the study by Cheng et al.⁷⁴ the maximum values of induced IFN-γ and IL2 were 4.7 × 10⁵ ng/ml and 5.2 × 10⁵ ng/ml, respectively. In the study by Binda et al.⁷⁵ the maximum values of induced IFN-γ and IL2 were 4.1 × 10⁵ ng/ml and 4.6 × 10⁵ ng/ml, respectively. We conducted four simulation rounds with different intervals to assess the impact of immunization schedules on immune responses. While innate immune responses remained consistent across schedules, adaptive immune responses varied. Memory B cells, activated B cells, and total TH cells peaked higher with a 1-week interval, while active TH cells peaked more with a 2-week interval. Cytokine production was schedule-dependent, with the 1-week interval schedule resulting in the highest peak values for IFN-γ, IL2, and TGF-β. Shorter intervals initiated the next round before cytokine levels fully declined, enhancing existing levels. Based on these findings, we opted to immunize mice with three injections at one-week intervals for animal experiments.

Delivery methods also play a significant role in vaccine efficacy. Subcutaneous vaccine administration has shown enhanced transgene expression and stronger antigen-specific CD8+ T cell-mediated immune responses compared to intramuscular injection.⁷⁶ Given the abundance of antigen-presenting cells in the epidermis and dermis, our study utilized subcutaneous injection with the MicronJet600 microneedle, as reported in previous studies,⁷⁷ to elicit superior immune responses and significant dose-saving effects after intradermal immunization.

In our animal model experiments, the candidate vaccine pVAX1-BD3-12P demonstrated robust immunogenic, therapeutic, and prophylactic effects. It is well known that protective responses are mediated by CD8 T cells in TC-1 tumor models,^78–80 and it has been demonstrated that CD4 T cells contribute to the effective activation of CD8 T cells and the development of long-term immune memory.^81,82 In turn, the novel multiepitope DNA vaccine developed in our study, encompassing multiple CTL and HTL epitopes, induced concurrent responses from CD4 and CD8 T cells. Our research reveals that our multi-epitope HPV vaccine not only achieves a good preventive effect but also demonstrates significant therapeutic efficacy. In our study, approximately 43% of mice treated with the vaccine experienced complete tumor regression, while the remaining mice showed inhibited tumor growth. Furthermore, in separate experiments, around 50% of vaccinated mice were still alive 60 days post-vaccination, compared to the control group mice euthanized within 35 days post-vaccination. Contrastingly, prior studies by Peng et al.^66,76 observed improved survival rates in tumor-bearing mice with codon-optimized DNA vaccines and needle-free delivery systems, respectively, yet complete tumor eradication was not achieved. Similarly, studies by Han et al.⁸³ found that their therapeutic vaccine alone showed a trend of inhibiting tumor growth but failed to completely eradicate tumors, whereas its effectiveness significantly improved when combined with Anti-PD-1/PD-L1. Additionally, studies by Matin Kayyal et al.⁸⁴ and Wang Qi et al.⁸⁵ reported significant inhibition of tumor growth with their multi-epitope vaccines, albeit without complete prevention or elimination of tumors. In summary, our research and previous DNA vaccine studies collectively suggest potential efficacy against HPV-related tumors, with variations in preventive effects, therapeutic effects, and survival times. Factors such as study design, animal models, and experimental conditions may influence specific immunogenicity and anti-tumor activity.

While our study provides valuable insights into the design of multi-epitope vaccines targeting HPV16, there are significant limitations that warrant attention. Despite in vitro simulations suggesting the safety of our vaccine, comprehensive in vivo and in vitro experimental validation is necessary to fully assess its safety profile. Additionally, experimental verification is currently lacking regarding the vaccine’s ability to activate TLR2, TLR3, and TLR4, as well as the impact of the β-defensin peptide on its immunogenicity. Furthermore, we need to consider conducting further experiments using gene knockout mice in the future to validate the precise role of CD4+ cells in our study. Filling these gaps is crucial for understanding the immunogenic mechanisms of the vaccine and its potential to induce robust immune responses. This will provide essential information for the future development of effective HPV vaccines.

Author contribution

L. Z., Y. T., X. C., L. S. (Lei Shi), X. Z., and Z. Y. made the contribution in acquisition of data, analysis and interpretation; L. Z. write the draft manuscript; Y. Y. and L. S. (Li Shi) design the study; L. S. (Li Shi) revised the manuscript. All authors read and approved the final manuscript.

Disclosure statement

We confirm that all of the listed authors actively participated in the present study and have read and approved the submitted manuscript. The authors do not have any potential conflicts of interest to declare.

Data availability statement

The authors confirm that the data supporting the findings of this study are available within the article or its supplementary materials.

Supplementary material

Supplemental data for this article can be accessed on the publisher’s website at https://doi.org/10.1080/21645515.2024.2352908

Additional information

Funding

This work was supported by Chinese Academy Medical Sciences Innovation Fund for Medical Sciences [2021-I2M-1-004], the Yunnan Provincial Science and Technology Department [202002AA100009], Yunnan Province Xingdian Talent Support Program (2024), the Yunnan Provincial Innovation Team of Therapeutic Neutralizing Antibody [202405AS350026] and NMPA Key Laboratory for Quality Control and Evaluation of Vaccines and Biological Products. The funders had no role in study design, data collection and analysis, decision to publish or preparation of the manuscript.