In silico/In vivo analysis of high-risk papillomavirus L1 and L2 conserved sequences for development of cross-subtype prophylactic vaccine

Human papillomavirus (HPV) is the most common sexually transmitted infection in the world and the main cause of cervical cancer. Nowadays, the virus-like particles (VLPs) based on L1 proteins have been considered as the best candidate for vaccine development against HPV infections. Two commercial HPV (Gardasil and Cervarix) are available. These HPV VLP vaccines induce genotype-limited protection. The major impediments such as economic barriers especially gaps in financing obstructed the optimal delivery of vaccines in developing countries. Thus, many efforts are underway to develop the next generation of vaccines against other types of high-risk HPV. In this study, we developed DNA constructs (based on L1 and L2 genes) that were potentially immunogenic and highly conserved among the high-risk HPV types. The framework of analysis include (1) B-cell epitope mapping, (2) T-cell epitope mapping (i.e., CD4+ and CD8+ T cells), (3) allergenicity assessment, (4) tap transport and proteasomal cleavage, (5) population coverage, (6) global and template-based docking, and (7) data collection, analysis, and design of the L1 and L2 DNA constructs. Our data indicated the 8-epitope candidates for helper T-cell and CTL in L1 and L2 sequences. For the L1 and L2 constructs, combination of these peptides in a single universal vaccine could involve all world population by the rate of 95.55% and 96.33%, respectively. In vitro studies showed high expression rates of multiepitope L1 (~57.86%) and L2 (~68.42%) DNA constructs in HEK-293T cells. Moreover, in vivo studies indicated that the combination of L1 and L2 DNA constructs without any adjuvant or delivery system induced effective immune responses, and protected mice against C3 tumor cells (the percentage of tumor-free mice: ~66.67%). Thus, the designed L1 and L2 DNA constructs would represent promising applications for HPV vaccine development.


Prediction of T-cell epitopes.
Since a linear form of T-cell epitopes are bound to MHCs, the interface between T-cells and ligands can be accurately modeled. In this study, we used three different algorithms (published motifs, ANN and quantitative matrix) for MHC-I and two algorithms for MHC-II (ANN and quantitative matrix).

Evaluation of L1 and L2 DNA expression in HEK-293T cells. In vitro DNA delivery of L1 and L2
into the eukaryotic cell line (HEK-293T) was performed by TurboFect as a transfection reagent. The levels of DNA expression were evaluated using fluorescence microscopy and flow cytometry at 48 h post-transfection. The data indicated that pEGFP-L1 and pEGFP-L2 can effectively penetrate into HEK-293T cells in vitro. The cellular uptake of the L1 and L2 genes into the HEK-293T cells was ~57.86% and ~68.42%, respectively. The delivery of pEGFP-N1 as a positive control was detected in approximately ~92.10% of HEK-293T cells (Fig. 3). Moreover, the spreading green regions were observed for L1 and L2 DNA delivery using TurboFect carrier by fluorescent microscopy in HEK-293T cells. On the other hand, western blot analysis indicated the successful expression of L1 and L2 proteins fused to GFP (i.e., L1-GFP and L2-GFP) using anti-GFP antibody. The data indicated the clear bands of ~52, ~50 and ~27 kDa for L1-GFP, L2-GFP and GFP, respectively using DAB substrate (Fig. 4).

Measurement of tumor growth.
To evaluate the prophylactic effects of the designed L1 and L2 DNA constructs, tumor growth and survival percentage were assessed in all groups for 60 days after challenging with C3 tumor cells. As shown in Fig. 5A, all test groups immunized with DNA constructs (G1, G2 & G3) demonstrated significantly lower tumor growth than that in control groups (PBS and empty vector, G4 & G5, p < 0.05). Our data showed progressive tumor growth in control groups on approximately 7-21 days (survival rate or tumor-free mice percentage: 0%). It was interesting that groups vaccinated with L1 DNA, L2 DNA and L1 + L2 DNA constructs similarly reduced the tumor growth (p > 0.05). As shown in Fig. 5B, group vaccinated with the mixture of L1 + L2 DNA constructs showed a higher survival rate (G3, ~66.67%) than L1 and L2 DNA constructs, alone (G1 & G2, ~33.33%).
Antibody assay. The levels of total immunoglobulin G (IgG), IgG2a and IgG2b in mice immunized with the mixture of L1 + L2 DNA constructs (G3) were significantly higher than other groups (p < 0.05, Fig. 6A,C,D).  Moreover, our data showed that the levels of IgG1 were similar in all groups vaccinated with DNA constructs (G1, G2 & G3, p > 0.05, Fig. 6B). There are no significant differences in the secretion of IgG2a and IgG2b isotypes between groups receiving the L1 and L2 DNA constructs, alone (G1 & G2, p > 0.05, Fig. 6C,D). No significant anti-(L1 + L2) antibody responses could be detected in the sera of control groups, thus, the seroreactivities were completely L1 + L2 antigen-specific responses in mice.
Granzyme B secretion. The secretion of Granzyme B in all test groups was significantly higher than the control groups (p < 0.05, Fig. 8). The group immunized with the L1 + L2 DNA construct (G3) produced significantly higher concentrations of Granzyme B than other groups (G1 & G2, p < 0.001). The level of Granzyme B in group receiving L1 DNA construct was similar to that in group receiving L2 DNA construct (p > 0.05).

Discussion
In recent years, development of bioinformatics tools applied in vaccine researches could potentially save time and resources. Indeed, the immunoinformatics tools help to identify antigenic domains for designing a multi-epitope vaccine. With sequence-based technology advancement, now we have enough information about the genomics and proteomics of different viruses 22 . Thus, using various bioinformatics tools, we can design peptide vaccines based on a neutralizing epitope. For example, in silico design of an epitope-based vaccine against human immunodeficiency virus 23,24 , coronavirus 25 , dengue virus 26 , and Saint Louis encephalitis virus 27 has already been reported.
While around 13 high-risk HPVs were recognized, current vaccines just protect humans from few types. An important limitation of the current vaccines is their narrow coverage. The accessibility of fully sequenced proteome from high-risk HPV strains provides a prospect for in silico screening of reliable peptide-based therapeutic vaccine candidates among billions of possible immunogenic peptides. In silico approaches are intended to reflect the possibilities for overcoming the above-mentioned difficulties in HPV multi-type vaccine. Gupta and coworkers designed prophylactic multiepitopic DNA vaccine using all the consensus epitopic sequences of HPVs L2 capsid protein. They also evaluated how engineering CpG motifs by bioinformatics tools could increase immunogenicity of DNA vaccines 28 . Hosseini et al. applied in silico analysis of L1 and L2 protein of HPV 11,16,18,31 and 45 types to identify universal peptide vaccine in order to protect against mentioned types 29 . In 2016, Singh et al. analyzed E1, E2, E6 and E7 proteins of high-risk HPV types to identify CD8 + T-cell epitopes. They suggested a pool of 14 peptides (9 to 43 amino acids) to provide the protection against high-risk HPV types 30 .  www.nature.com/scientificreports www.nature.com/scientificreports/ Panahi and colleagues used a two-step method (consist of molecular docking and sequence-based approach) to determine immunogenic epitopes for induction of immune system against the oncoproteins of HPV 16, 18, 31 and 45 types 31 34 .
In this research, we designed a framework for the comprehensive analysis of L1 and L2 conserved regions of high-risk HPV types containing both MHC-I and MHC-II epitopes. The framework begins with conservancy analysis of all 13 high-risk HPV strains following with (1) B-cell epitope mapping, (2) T-cell epitope mapping (CD4 + and CD8 + ), (3) allergenicity assessment, (4) tap transport and proteasomal cleavage, (5) population coverage, (6) global and template-based docking and (7) data collection, analysis, and design of the L1 and L2 DNA constructs. For experimental analysis, the final L1 or L2 DNA constructs were cloned into mammalian expression vector with green fluorescent tag (pEGFP vector) and their expression was evaluated in the eukaryotic cells using flow cytometry, fluorescent microscopy and western blotting. Moreover, the L1/L2-specific antibody and T-cell immune responses induced by L1 and L2 DNA constructs were assessed in mouse tumor model.
At first, L1 and L2 sequences obtained from high-risk HPV types were aligned using MUSCLE algorithms. Conservancy analysis showed that five regions of HPV16,18 L1 protein (8-22, 95-132, 307-342, 398-425 and 449-473) and four regions of HPV16,18 L2 protein (11-40, 54-76, 96-120 and 278-305) were more conserved among other subtypes and could be analyzed as an immunoinformatics input. In B-cell epitope prediction, L1 [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22] , L1 408-421 , L1 404-417 , L2 22-35 , L2 100-113 , L2 94-107 and L2 57-70 had the highest epitope prediction scores. Unfortunately, a reliable method for prediction of B-cell epitope has not been revealed up to now and the sensitivity and specificity of existing methods were very low (the specificity and sensitivity of this method were 0.57 and 0.58, respectively). In the case of T-cell epitope prediction, in silico analysis has been significantly improved, thus, the results are more reliable. In this study, for MHC-I epitopes, L1 12-21 (YLPPVPVSKV-type16 and YLPPPSVARV-type18), L1 460-470 (DQFPLGRKFLL-type16), L1 461-471 (DQYPLGRKFLV-type18), L2 11-20 (KRASATQLYK-type16 and KRASVTDLYK-type18), L2 280-291 (DPDFLDIVALHR-type16) and L2 273-284 (DSDFMDIIRLHR-type18) epitopes had the highest binding affinity scores. In addition, above-mentioned epitopes had the highest T-cell epitope prediction scores which were obtained from proteasomal cleavage and tap transport analysis. High degree of conservancy was observed between subtypes for these epitopes (  [54][55][56][57][58][59][60][61][62][63][64][65][66][67][68][69] (FFGGLGIGTGSGTGGR-type16) epitopes had the highest binding affinity scores. Among them, L2 54-69 had the greatest degree of conservancy (high similarity with all of the high-risk HPV types). One of the remarkable points is that L1 [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22] and L2 57-70 epitopes are the same (or overlapping with little difference (among B-cell and MHC-II selected epitopes. Due to a limitation of MHC-peptide binding prediction such as the gap between the peptides that are predicted to bind to MHC and those that experimentally bind 35   www.nature.com/scientificreports www.nature.com/scientificreports/ been employed to address this problem and raise the accuracy of MHC-peptide prediction. In the current study, template-based docking and also global docking were performed on the selected peptides to determine which peptide would get into the groove of MHC with the highest modeling scores. For MHC-I epitope, L1 12-21 , L1 460-470 and L2 280-291 sequences had the highest interaction similarity and cluster density scores. For MHC-II epitopes, L1 95-111 , L1 417-431 , L2 100-118 and L2 281-297 sequences had the highest docking scores. In this study, MHC-I-peptide docking scores confirmed MHC-I-peptide binding affinity scores because the same epitopes had the highest scores in both methods but in MHC-II molecular docking, the results were slightly different. One of the reasons is the significant conformational changes during the process due to the longer epitope length. As a general rule: the longer the length of the query peptide, the more torsions and conformational flexibilities 36 . Herein, due to longer peptide sequences, docking results in MHC-II were less accurate than MHC-I. For example, average similarity score in MHC-I was variable (171.8-259.7), but in MHC-II was 115.4-136. After the completion of the analysis and according to all of the above-mentioned parameters, two separate constructs were designed. In addition, accumulative population coverage of helper T-cell and CTL epitopes for the designed constructs   www.nature.com/scientificreports www.nature.com/scientificreports/ were estimated. For the L1 and L2 constructs, the combination of 8 epitope candidates for helper T-cell and CTL in a single universal vaccine could involve all world population by the rate of 95.55% and 96.33%, respectively (Fig. 3). In previous studies, YLPPVPVSKV (HPV16 L1) 37 and KRASVTDLYK (HPV18 L2) 21 have been reported as potentially immunogenic epitopes. The ability of in vitro expression of the designed L1 and L2 DNA constructs was determined in HEK-293T cells using flow cytometry and western blot analysis. The transfection efficiency of the L1 and L2 DNA constructs was ~57.86% and ~68.42%, respectively indicating their high potency for delivery into the eukaryotic cells. As known, the use of a polytope DNA vaccine containing multiple T-cell and B-cell epitopes is an attractive strategy for developing a therapeutic and prophylactic vaccine against HPV infections. After in vitro assay, immunological experiments were performed in mice to determine the efficiency of the designed L1 and L2 DNA constructs without the use of adjuvant or delivery system for vaccine development. Similarly, some studies used the pcDNA vector harboring the gene of interest for immunization without any adjuvant 38,39 . Our data indicated that the groups immunized with L1, L2 and L1 + L2 DNA constructs increased antibody and T-cell responses as compared to control groups. Furthermore, the (L1 + L2)-specific immunity in mice receiving the mixture of L1 + L2 DNA constructs (G3) resulted in higher secretion of total IgG, IgG2a, IgG2b, IFN-γ, IL-5 and IL-10 cytokines as well as Granzyme B than other groups. The higher levels of IgG2a and IgG2b as well as IFN-gamma (as a Th1 cytokine) in this group drive T-cell responses toward Th1-type immunity. The studies showed that immunoglobulin G1 (IgG1) is related to a Th2-type response, while a Th1 response is associated with the induction of IgG2a and IgG2b in mice 40 . Regarding to our observations in protective studies, this regimen (L1 + L2 DNA construct: G3) could confer further protection against C3 tumor-challenged mice (survival rate: ~66.67%) depending on stimulation of CD4 + T cell-dominated Th1 responses as well as Granzyme B secretion (indicating CTL activity) as compared to the L1 or L2 DNA constructs, alone (survival rate: ~33.33%). These data showed high potency of the combined L1 + L2 DNA constructs versus each DNA construct alone as a prophylactic HPV vaccine. Taken together, immunoinformatics approaches have been emerged as a critical field for accelerating immunological researches. Yet, the immunoinformatics techniques applied to T-cells have more advancement than those dealing with B-cells 30 . Moreover, recently, due to the limited options for choosing an adjuvant in clinical trials, bioinformatics analyses have been developed to predict the best adjuvant. In this way, in silico studies help researchers saving time and resources, and also can guide the experimental work with higher probabilities of finding the desired solutions and with fewer trial and error repeats of assays. The accessibility of HPV genomic sequences and functional characterization of the genes involved in the virulence has significantly improved our understanding of the molecular foundation for the pathogenesis of HPV and offered a wealth of data that can be used to design new plans for vaccine design. Nowadays, powerful immune system simulators have been developed using bioinformatics tools which predict artificial immunity provided by the vaccine. These approaches could predict the best adjuvant for using in human vaccine studies. There is a multi-scale computational infrastructure approach which can stimulate the dynamics of the immune response induced by several vaccination formulations and predict optimal combination in terms of adjuvant type, dosage and timing. NetLogo is an agent-based modeling of the immune system running different simulations with different parameter settings. It also can interact with different modeling strategies including the investigation of pathogen growth, life cycle modeling environment for simulation complex phenomena [41][42][43] . Therefore, using these methods can increase efficiency and reduce costs in vaccine studies. In this study, for the first time, comprehensively integrated methods (using sequence-based tools in combination with flexible peptide-protein docking) were used to design highly immunogenic and protective vaccine candidates which were able to boost both humoral and cellular    Table 12. MHC-II -peptide docking scores of selected helper T-cell epitopes. *higher rate shows better quality of peptide-MHC interactions. www.nature.com/scientificreports www.nature.com/scientificreports/ immune responses against all high-risk HPV types. In addition, in vivo analysis demonstrated high potency of the designed L1 and L2 constructs as combined in DNA-based vaccines without the use of adjuvant or delivery system. However, we will improve the efficiency of these DNA-based vaccines using a delivery system and also will compare their efficacy with the designed peptide-based vaccines along with adjuvants in near Future.      Table 15. Physicochemical properties of L1 and L2 DNA vaccine constructs. *higher rate shows high degree of peptide antigenicity. **higher rate shows high degree of peptide solubility. Protein alignments and conservancy analysis. To determine conserved epitopes between different subtypes, L1 and L2 sequence datasets were first aligned using SnapGene software 4.2.2 (From GSL Biotech; available at snapgene.com). After protein alignments analysis using muscle algorithms, the conserved epitopes of each protein were selected for immune-bioinformatics analysis such as B-and T-cell epitope prediction. Also, to calculate the degree of variability and conservancy of each epitope, IEDB epitope conservancy tools (http://tools.immuneepitope.org/tools/conservancy/) were used.

MHC-I MHC-II
Linear B-cell epitope prediction. A successful vaccine must elicit a strong T-cell and B-cell immune response, but above all, provide protection against the disease being targeted. Therefore, it is essential to show that constructed immunogens are able to induce protective cellular and humoral immunity. Since the antibodies are induced against linear B-cell epitopes, it would be very difficult to synthesize long peptides with the native protein conformation resembling for the induction of protective antibodies. However, optimal peptide-based vaccines should be presented in a desired secondary structure of peptides in order to induce a specific humoral response 41,42 . For the B-cell epitope prediction of conserved regions in L1 and L2 proteins, BepiPred-2.0 server (http://www.cbs.  www.nature.com/scientificreports www.nature.com/scientificreports/ dtu.dk/services/BepiPred-2.0/) was employed. In this study, epitope threshold value was set as 0.5 (the specificity and sensitivity of this method are 0.57 and 0.58, respectively) 41 .
T-cell epitope prediction. MHC-I epitope prediction: The initial step on applying bioinformatics to vaccine researches is to assess potentially immunoprotective epitopes. T-cell epitopes presented by MHC molecules are typically in a linear form containing 12 to 20 amino acids. This fact facilitates accurate modeling for the interaction of ligands and T-cells 44 . Thus, the most selective step in the presentation of antigenic peptide to T-cell receptor (TCR) is the binding of the MHC molecule 45 . In this study, we tried to use three different algorithms including Artificial Neural Networks (NetMHCpan 4.0 server 43 (http://www.cbs.dtu.dk/services/NetMHCpan/), Quantitative matrix (Propred I 43 (http://crdd.osdd.net/raghava/propred1/) and Published motifs (syfpeithi server 46 (http://www.syfpeithi.de) to predict high-potential T-cell epitopes. For NetMHCpan, percentile rank was set at 0.5% for strong binders and 2% for weak binders and for Propred I threshold was set at 4%.
Prediction of MHC-I peptide presentation pathway. Investigating the Tap transport and proteasomal cleavage as well as affinity prediction of binding is essential in MHC-I presentation pathway. In this study, we used NetCTL 1.2 server combined with Tap transport/proteasomal cleavage tools (http://www.cbs.dtu.dk/services/NetCTL/) to access the prediction of antigen processing through the MHC class I antigen presentation pathway. In this method, parameters of weight on the C-terminal cleavage, Tap transport efficiency, and epitope identification were set to default (0.15, 0.05 and 0.75, respectively) 49 .
Population coverage. Since the response to T-cell epitopes is restricted by MHCs, the selection of epitopes with multiple HLA-binding increases population coverage in defined geographical regions where the peptide-based vaccine might be employed. The coverage rate of population for each epitope was computationally validated using the IEDB population coverage tool 50 (/population/iedb_input). In this study, individual epitope and its binding to HLA alleles were analyzed, and different geographic areas were also selected.
Allergenicity and cross-reactivity assessment. Since proteins are very important in inducing allergenic reactions, the prediction of potential allergenicity is an important item in the safety assessment especially in the field of genetically modified foods, therapeutics, bio-pharmaceuticals etc. 51 . The food and agriculture organization (FAO) and world health organization (WHO) protocol includes three terms to evaluate the allergenicity of proteins which are defined as following: the term sensitivity refers to correctly predicted allergens (%), whereas www.nature.com/scientificreports www.nature.com/scientificreports/ specificity refers to correctly predicted non-allergens (%), and also accuracy refers to the proportion of correctly predicted proteins 19 . The allergenicity of the epitopes was analyzed by the PA 3 P (http://lpa.saogabriel.unipampa. edu.br:8080/pa3p/pa3p/pa3p.jsp) using Allergen online (8aa and 80 wordmatch) and AFDS-motif algorithms based on amino acid composition. The specificity of these methods is 95.43% (8aa), 92.88% (80aa) and 88.1% (ADFS) 52 . To assess cross-reactivity between peptide and human proteome, top-ranked epitope were analyzed by peptide matching program (https://research.bioinformatics.udel.edu/peptidematch/index.jsp) 53 .
Peptide-protein flexible Docking. Computational docking methods have been known as an important tool for drug design 54 . With the rapid development of peptide therapeutics in rational drug design, the use of new techniques such as protein-peptide docking is inevitable. In this study, two different algorithms (template-based docking and global docking) were performed by GalexyPepDock server 55 (http://galaxy.seoklab.org/cgi-bin/ submit.cgi?type=PEPDOCK) and CABS Dock server 56 (http://biocomp.chem.uw.edu.pl/CABSdock). To estimate the formation of MHC-peptide complex, the GalaxyPepDock server effectively models the structural 3D peptide-protein complexes from input peptide and protein sequences using the structure database and energy-based optimization (Template-based Docking). CABS-Dock server performs Global docking procedure which at first explicit fully flexible docking simulation and then clustering-based scoring. Receptor flexibility was limited by default to small backbone fluctuation but could be increased to include selected receptor fragments 56,57 . This study presented an example of MHC-peptide docking performed by each individual epitope and available PDB file (Table 13) of HLA alleles, separately.
Physicochemical properties of the designed L1 and L2 constructs. Based on L1 and L2 top-ranked epitopes, two different constructs were designed. The physicochemical properties of top-ranked epitopes such as solubility, molecular weight, estimated half-time, instability index and antigenicity were determined by ProtParam (https:// web.expasy.org/protparam/) tools 58 , VaxiJen 59 (http://www.ddg-pharmfac.net/vaxijen/VaxiJen/VaxiJen.html) and Protein-Sol (https://protein-sol.manchester.ac.uk/) server 60 . www.nature.com/scientificreports www.nature.com/scientificreports/ Experimental studies Construction of the recombinant plasmids. After bioinformatics analysis, the selected peptides were assembled in two separated constructs (Fig. 2). The pUC57-L1 and pUC57-L2 constructs were synthesized by Biomatik Company. For in vitro experiments, the pUC57-L1 and pUC57-L2 vectors were digested by XhoI/HindIII, and the L1 and L2 genes were subcloned into XhoI/HindIII sites of pEGFP-N1 vector, individually (i.e., pEGFP-L1 and pEGFP-L2). All the recombinant vectors were transformed into Escherichia coli (E.  www.nature.com/scientificreports www.nature.com/scientificreports/ coli) DH5α strain. After extraction of plasmids from single colonies using Mini-Kit (Qiagen), the presence of inserted L1 and L2 fragments was confirmed by digestion with restriction enzymes and sequencing. For in vivo immunological assessment, the pUC57-L1 and pUC57-L2 vectors were digested by BamHI/HindIII and the L1 and L2 genes were subcloned into BamHI/HindIII sites of pcDNA3.1 (-) vector containing cytomegalovirus early promoter and enhancer sequence, individually (i.e., pcDNA-L1 and pcDNA-L2). Indeed, we used the pcDNA vector harboring CpG motif for in vivo studies. As a final point, the recombinant DNA vectors harboring L1 and L2 genes were purified by an endotoxin-free plasmid Extra EF kit (Macherey Nagel, Germany). The concentration and purity of the recombinant L1 and L2 DNA constructs were determined by NanoDrop spectrophotometry 61 . In vitro expression of L1 and L2 DNA constructs in HEK-293T cells. Human embryonic kidney cells (HEK-293T) were cultured in RPMI supplemented with 10% fetal bovine serum (FBS) at 37 °C and 5% CO 2 atmosphere. After some passages, the cells were seeded in a 12-well plate. The optimal cell confluency for effective transfection was considered 70-80%. For the generation of TurboFect-plasmid DNA complex, 10 μl of TurboFect (Thermo Scientific) and 2 μg of each plasmid (pEGFP-L1, pEGFP-L2 and pEGFP-N1 as a positive control) were mixed and incubated for 15 min at room temperature. Then, the complex was added to each well in serum-free media. In addition, the non-transfected HEK-293T cells were used as negative control. After six hours, the media was replaced with the completed RPMI medium. Finally, the cells were harvested, washed and resuspended in PBS buffer, to analyze the expression of L1 and L2 DNA constructs using flow cytometry, fluorescent microscopy and western blotting at 48 hr after transfection 61 .
Mice immunization. Five groups of six female C57BL/6 mice (obtained from the breeding stocks maintained at Pasteur Institute of Iran; MHC haplotype B/H-2Kb/H-2Db) were immunized on days 0, 14, and 28 (i.e., three times with a 2-week interval) with 50 µg of each plasmid DNA (pcDNA-L1 or pcDNA-L2: G1 or G2) or their combination (pcDNA-L1+ pcDNA-L2: G3) at the right footpad as shown in Table 16. The control groups (G4 and G5) received pcDNA3.1 and PBS, respectively. All mice were maintained under specific pathogen-free conditions 62 . Moreover, all of the animal experimental procedures were approved by Animal Care and Use Committee of Pasteur Institute of Iran and carried out according to the Animal Experimentation Regulations of Pasteur Institute of Iran (national guideline) for scientific purposes (code: 976).

Monitoring tumor growth.
For in vivo protection assay, vaccinated mice were subcutaneously challenged in the right flank with C3 tumor cells (1 × 10 5 cells), two weeks after the last injection. The C3 tumor cells contain whole HPV16 genome, and the presence of L1 and L2 genes was confirmed in the previous studies 63 . Tumor growth and the percentage of tumor-free mice were monitored twice a week by palpation for 60 days post-challenge. At each time, tumor volume was calculated by this formula: V = (a 2 b)/2 (a = the smallest diameter and b = the biggest diameter) 62 .
Antibody assay secreted from B-cells. Two weeks after the last injection, serum samples were collected from each group. The levels of goat anti-mouse immunoglobulin G1 (IgG1), IgG2a, IgG2b and total IgG antibodies (diluted 1:10,000 in 1% BSA/PBS-Tween, Sigma) secreted from B-cells were measured in the pooled sera of each group by indirect ELISA. The coated antigens were the mixture of L1 and L2 synthetic peptides (5 μg/mL). Moreover, mice sera were diluted 1:100 in 1% BSA/PBS-Tween 64 .
Cytokine assay secreted from T-cells. Three mice from each group were sacrificed and the spleens were removed. The red blood cell-depleted pooled splenocytes (2 × 10 6 cells/ml) were cultured in 48-well plates for 72 h in the presence of 5 μg/mL of L1 + L2 peptides, RPMI 5% as negative control and 5 μg/mL of concanavalin A (ConA) as positive control in complete RPMI culture medium. The supernatants were harvested to assess the secretion of IFN-γ, IL-5 and IL-10 from T-cells using the sandwich-based ELISA method (R&D Systems) according to the manufacturer's instructions. All data were represented as mean ± SD for each sample 65 . Granzyme B assay (in vitro CTL activity). To measure Granzyme B (GrB) by ELISA, the P815 target cells (T) were seeded into 96-well plates (2 × 10 4 cells/well) incubated with the mixture of L1 and L2 peptides (~30 μg/mL) for 24 h. Then, the prepared splenocytes (Effector cells: E, before section) were counted and added to the target cells at E: T ratio of 100: 1 in complete RPMI culture medium for 6 h incubation. Finally, the supernatants were harvested to measure the concentration of GrB by ELISA (eBioscience kit) according to the manufacturer's instruction 64 .  Table 16. Immunization program for in vivo analysis.