key: cord-340781-z348xbn0
authors: Namvar, Ali; Bolhassani, Azam; Javadi, Gholamreza; Noormohammadi, Zahra
title: In silico/In vivo analysis of high-risk papillomavirus L1 and L2 conserved sequences for development of cross-subtype prophylactic vaccine
date: 2019-10-23
journal: Sci Rep
DOI: 10.1038/s41598-019-51679-8
sha: 
doc_id: 340781
cord_uid: z348xbn0

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.

18 31 www.nature.com/scientificreports www.nature.com/scientificreports/ charge and secondary structure. At first step, the conserved region sequences were analyzed by BepiPred-2 server to predict potential B-cell epitopes (Table 3 ). In L1 protein, L1 [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] (EATVYLPPVPVSKVV-type16), L1 408-421 (PPPGGTLEDTYRFV-type16) and L1 404-417 (NFGVPPPPTTSLVD-type 18) epitopes had the best B cell epitope identification scores. For L2 protein, L2 [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] (KQSGTCPPDVVPKV-type18), L2 100-113 (PSDPSIVSLVEETS-type16), L2 94-107 (EPVGPTDPSIVTLI-type18) and L2 [57] [58] [59] [60] [61] [62] [63] [64] [65] [66] [67] [68] [69] [70] (GLGIGTGSGTGGRT-type16) epitopes showed the highest epitope identification score between their own protein sequences.

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).

Prediction of MHC-I. At first step, the L1 and L2 conserved regions were analyzed to find the most immunodominant peptides using NetMHCpan 4.0, syfpeithi and ProPred I. In each protein, peptides with the highest binding affinity scores were determined as high-potential CTL epitope candidates (Tables 4 and 5 ). The analysis showed that L1 [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] (YLPPVPVSKV-type 16 and YLPPPSVARV-type 18), L1 460-470 (DQFPLGRKFLL-type 16) , L1 461-471 (DQYPLGRKFLV-type 18), L2 [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] (KRASATQLYK-type 16 and KRASVTDLYK-type 18), L2 280-291 (DPDFLDIVALHR-type 16) and L2 273-284 (DSDFMDIIRLHR-type 18) epitopes had the highest binding affinity among their own protein sequences. In general, the results of three different algorithms confirmed each other. Conservancy and allergenicity analyses were done on the selected epitopes. The sequence of all the epitopes were well conserved among high-risk HPV types and none of them were allergens (Tables 4 and 5 ). In addition, there was no cross-reactivity between peptide and human proteome.

In this study, we used NetMHCIIpan and Propred servers for MHC-II epitope identification analysis (Table 6 ). Since a suitable T-cell epitope should be predicted to bind to different HLA alleles, epitopes with the maximum number of binding HLA-DR alleles were selected as high-potential helper T-cell epitope candidates. Among predicted epitopes, L1 [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] (EATVYLPPVPVSKVV-type 16) , L1 95-111 (TQRLVWACVGVEVGRGQ-type 16 and TQRLVWACAGVEIGRGQ-type 18), L1 416-430 (DTYRFVTSQAIACQK-type 16) , L1 417-431 (DTYRFVQSVAITCQK-type 18), L2 100-118 (DPSIVTLIEDSSVVTSGAP-type 16) , L2 281-297 (PDFLDIVALHRPALTSR-type 16) and L2 274-290 (SDFMDIIRLHRPALTSR-type 18) had the highest scores of binding affinity. Also, the sequence of all the epitopes were well conserved among high-risk HPV types and none of them were allergen (Tables 6 and 7) . Also, there was no cross-reactivity between peptide and human proteome. Population coverage analysis. HLA distribution varies among the diverse geographic regions around the world. Thus, while designing an effective vaccine, population coverage must be taken into consideration to cover the maximum possible populations. In this study, population coverage was estimated separately for each putative epitope in 16 specified geographic regions of the world (Tables 9 and 10) . For CTL epitopes, the highest population coverage of world's population was calculated for L1 12 a clear band of ~765 bp and ~700 bp on agarose gel for L1 and L2, respectively (data not shown). The recombinant endotoxin-free plasmids (i.e., pcDNA-L1 and pcDNA-L2) had a concentration range between 1.5 and 3.5 mg/mL.

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 ).

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.

Cytokine assay. The results of cytokine assay in each group showed that the levels of (L1 + L2)-specific IFN-γ, IL-10 and IL-5 secretions in groups immunized with L1 (G1), L2 (G2) and L1 + L2 (G3) DNA constructs were significantly higher than control groups (p < 0.05, Fig. 7 www.nature.com/scientificreports www.nature.com/scientificreports/ Fig. 7) . Furthermore, our data indicated that the ratios of IFN-γ/IL-10 and IFN-γ/IL-5 were higher in all test groups as compared to control groups; therefore, they could trigger Th1 immune response.

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).

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.

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%.

MHC-II epitope prediction: For MHC class II, NetMHCIIpan 3.2 server 47 (http://www.cbs.dtu.dk/services/ NetMHCIIpan/) and ProPred 48 (http://crdd.osdd.net/raghava/propred/) were employed to predict potential interaction of helper T-cell epitope peptides and MHC-II. In this case, the threshold for strong and weak binders was set at 2% and 10%, respectively.

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 .

Western blot analysis. HEK-293T cells were scraped from their plates and washed with PBS1X. After washing steps, the cells were lysed in whole-cell lysis buffer (10% glycerol, 1 nM DTT, 2 mM natrium fluoride, 0.2% Triton X-100, 0.5 EDTA in PBS pH = 7.4). The extracted protein samples (L1-GFP, L2-GFP and GFP) were separated by SDS-PAGE in 12.5% (w/v) polyacrylamide gel and transferred to nitrocellulose membrane (Millipore). The membrane was equilibrated with TBST (Tris-buffered saline Tween-20) solution containing 2.5% BSA (Bovine albumin serum) overnight. The anti-GFP polyclonal antibody (1:5000 v/v; Acris antibodies GmbH) was used to recognize the expressed proteins under standard procedures. The immunoreactive protein bands were visualized by detection of peroxidase activity using a substrate named as 3, 3′-diaminobenzidine (DAB, Sigma) 61 .

Peptide constructs synthesis. For immunological assay (i.e., secretion of antibody, cytokine and Granzyme B), two peptide constructs (L1 and L2 peptides, Fig. 2) were synthesized by BioMatik Co. with more than 85% purity.

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).

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.

Worldwide burden of cancer attributable to HPV by site, country and HPV type

Estimate of the global burden of cervical adenocarcinoma and potential impact of prophylactic human papillomavirus vaccination

HPV vaccination: the promise & problems

Pros, cons, and ethics of HPV vaccine in teens-Why such controversy?

Virus-like particles for the prevention of human papillomavirus-associated malignancies

Therapeutic human papillomavirus vaccines: current clinical trials and future directions

Efficacy of fewer than three doses of an HPV-16/18 AS04-adjuvanted vaccine: combined analysis of data from the Costa Rica Vaccine and PATRICIA Trials

VLPs displaying a single L2 epitope induce broadly crossneutralizing antibodies against human papillomavirus

The papillomavirus virion: a machine built to hide molecular Achilles' heels

Concatenated multitype L2 fusion proteins as candidate prophylactic pan-human papillomavirus vaccines

Optimization of multimeric human papillomavirus L2 vaccines

Potent anti-HPV immune responses induced by tandem repeats of the HPV16 L2 (20-38) peptide displayed on bacterial thioredoxin

Chimeric L1-L2 virus-like particles as potential broad-spectrum human papillomavirus vaccines

Immunogenic assessment of plant-produced human papillomavirus type 16 L1/L2 chimaeras

Protection against heterologous human papillomavirus challenge by a synthetic lipopeptide vaccine containing a broadly cross-neutralizing epitope of L2

Generation and characterization of a preventive and therapeutic HPV DNA vaccine

Control of HPV-associated tumors by innovative therapeutic HPV DNA vaccine in the absence of CD4+ T cells

Cancer Immunoinformatics: A Promising Era in the Development of Peptide Vaccines for Human Papillomavirus-induced Cervical Cancer

Computational Perspective on the Current State of the Methods and New Challenges in Cancer Drug Discovery

Designing of CD8(+) and CD8(+)-overlapped CD4(+) epitope vaccine by targeting late and early proteins of human papillomavirus

An overview of bioinformatics tools for epitope prediction: Implications on vaccine development

Prediction of cross-clade HIV-1 T-cell epitopes using immunoinformatics analysis

Immune responses induced in HHD mice by multiepitope HIV vaccine based on cryptic epitope modification

Design of an epitope-based peptide vaccine against spike protein of human coronavirus: an in silico approach

A computational approach for identification of epitopes in dengue virus envelope protein: a step towards designing a universal dengue vaccine targeting endemic regions

A computational assay to design an epitope-based Peptide vaccine against saint louis encephalitis virus

In silico DNA vaccine designing against human papillomavirus (HPV) causing cervical cancer

Silico Analysis of L1/L2 Sequences of Human Papillomaviruses: Implication for Universal Vaccine Design

Sequence-based approach for rapid identification of cross-clade CD8+ T-cell vaccine candidates from all high-risk HPV strains. 3 Biotech 6

Identification and in silico analysis of functional SNPs of human TAGAP protein: A comprehensive study

Designing of Epitope-Focused Vaccine by Targeting E6 and E7 Conserved Protein Sequences: An Immuno-Informatics Approach in Human Papillomavirus 58 Isolates

Identification of a ras oncogene peptide that contains both CD4(+) and CD8(+) T cell epitopes in a nested configuration and elicits both T cell subset responses by peptide or DNA immunization

Prediction of MHC-peptide binding: a systematic and comprehensive overview

Modeling of protein-peptide interactions using the CABS-dock web server for binding site search and flexible docking

Interaction of the host immune system with tumor cells in human papillomavirus associated diseases

Multivalent human papillomavirus l1 DNA vaccination utilizing electroporation

Enhancing antitumor immunogenicity of HPV16-E7 DNA vaccine by fusing DNA encoding E7-antigenic peptide to DNA encoding capsid protein L1 of Bovine papillomavirus

Th1-directing adjuvants increase the immunogenicity of oligosaccharide-protein conjugate vaccines related to Streptococcus pneumoniae type 3

Combining agent based-models and virtual screening techniques to predict the best citrus-derived vaccine adjuvants against human papilloma virus

A computational model to predict the immune system activation by citrus-derived vaccine adjuvants

A growth model of human papillomavirus type 16 designed from cellular automata and agent-based models

T-cell epitope mapping

Immunodominance in major histocompatibility complex class I-restricted T lymphocyte responses

SYFPEITHI: database for MHC ligands and peptide motifs

Improved methods for predicting peptide binding affinity to MHC class II molecules

ProPred: prediction of HLA-DR binding sites

Large-scale validation of methods for cytotoxic T-lymphocyte epitope prediction

Predicting population coverage of T-cell epitope-based diagnostics and vaccines

AlgPred: prediction of allergenic proteins and mapping of IgE epitopes

Prediction of protein allergenicity based on signal-processing bioinformatics approach

A fast Peptide Match service for UniProt Knowledgebase

computer-enabled peptide drug design: principles, methods, applications and future directions

GalaxyPepDock: a protein-peptide docking tool based on interaction similarity and energy optimization

CABS-dock web server for the flexible docking of peptides to proteins without prior knowledge of the binding site

Immune Responses and Anti-Tumor Potential of an HPV16 E6E7 Multi-Epitope Vaccine

The Proteomics Protocols Handbook

VaxiJen: a server for prediction of protective antigens, tumour antigens and subunit vaccines

Protein-Sol: a web tool for predicting protein solubility from sequence

HPV16 L2 improves HPV16 L1 gene delivery as an important approach for vaccine design against cervical cancer

Whole recombinant Pichia pastoris expressing HPV16 L1 antigen is superior in inducing protection against tumor growth as compared to killed transgenic Leishmania

Immunogenicity of an HPV-16 L2 DNA vaccine

Small heat shock protein 27: An effective adjuvant for enhancement of HIV-1 Nef antigen-specific immunity

Recombinant Leishmania tarentolae encoding the HPV type 16 E7 gene in tumor mice model

A.N. and A.B. conceptualized the work. A.N. performed the experiments. A.N., A.B., G.J. and Z.N. analyzed the data. A.N. wrote the manuscript. A.B. edited the manuscript. All authors approved the final version of the paper.

The authors declare no competing interests.

Correspondence and requests for materials should be addressed to A.B.Reprints and permissions information is available at www.nature.com/reprints.Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.