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

Abstract

The aim of this study was to design a multiepitope universal vaccine for major human papillomavirus (HPV) structural proteins, L1/L2, by bioinformatics models. For this purpose, we predicted the most probable immunogenic epitopes of L1 and L2 from common high-risk HPV 16, 18, 31, and 45 beside high prevalent type 6 and 11 based on BCPREDS defaulted model, while solvent accessibility of structure was extrapolated. The three-dimensional molecular model of L1 protein was constructed by Swiss Model server, whereas sequence alignment provided model for prediction of L2 protein epitopes. After that, N-glycosylation sites were excluded from estimated epitope regions. Then, by other bioinformatics analyses, 20 epitopes were selected and fused in tandem repeats, reverse translated, and codon optimized to relevant sequence. The final protein parameters such as antigenicity were analyzed by protean program. Evaluation of new recombinant protein sequence indicated a molecular weight of 41.8 kDa with 400 amino acids beside positive charge. The computed isoelectric point (pI) value indicated the acidic nature of final product. The aliphatic index showed low thermal stability of this construct and the Grand Average Hydropathicity value was negative (−0.494). Analyzed plot showed that major parts of new protein construct had hydrophilic property, thus harboring antigenic potency. After all, sequence of final construct reverse translated to DNA and this codon-optimized sequence showed Codon Adaptation Index (CAI) of >0.8 for expression in Escherichia coli. Finally, this sequence ligated into pET28a bacterial expression vector. The new recombinant proteins harboring 20 B cell epitope seem to be suitable antigens based on computational methods as a universal vaccine candidate for HPVs.

Introduction

Human papillomavirus (HPV), as an epitheliotropic DNA virus, is a nonenveloped virion responsible for spectrum of benign (such as warts and papillomas) to malignant mucosal lesion in humans (4,18). More than hundreds of papillomavirus types have been distinguished based on L1 protein sequence in mammalian host, and among them, the small subset has medical importance in humans (12). Near 40 types of HPVs can establish a persistent infection inside the genital tract, but only a few ultimately lead to cervical cancer, one major cause of women death with around 530,000 new cases per year (1). The clinically important types are categorized into low-risk or high-risk groups by considering the relationship with malignancy. While infection with type 6 and 11 as low-risk HPVs causes less serious outcome such as genital warts, high-risk HPVs, including types 16, 18, 31, and 45, are responsible for >90% of cervical cancer worldwide (8,12,20).

The papillomavirus genome is ∼7.9–8 kb long and contains three main regions: a noncoding region (NCR), termed also as long-control region (LCR) or upstream regulatory region (URR), an early (E) region of open reading frames (ORFs), and a late (L) region of ORFs. The circular papillomavirus genome encodes roughly eight ORFs from a single DNA strand (18), which is responsible for viral replication, maintenance, and cell transformation. These gene products have been classified into six regulatory (E1-E6) and two structural proteins (L1 and L2) (3). An extensive body of researches has demonstrated that E1 and E2 are involved in viral transcription and replication; E5, E6, and E7 act as oncoproteins, while L1 and L2 contribute to the formation of the capsid (new). The viral icosahedral capsid is built up from 72 identical capsomers that comprised five L1 and one L2 molecule. The exposed L1 and L2 proteins have been identified as competent targets for prophylactic vaccination as they are expressed during early infection (13). The L1 is the major capsid protein with molecular weight of 55 kDa, being particularly abundant and representing 80% of the capsid proteins (ratio L1/L2 10:1). L2 protein is about 500 amino acids long with molecular weight of 64–78 kDa and representing ∼20% of the capsid proteins, thus being poorly represented on the surface of HPV virions (18).

Even providing an effective vaccine did not control the rate of cervical cancer in population of developing countries. Different kinds of vaccination strategies have been developed for HPV; among them, peptide-based and virus-like particle vaccines have considerable value (13). The FDA-approved HPV vaccines designated as Gardasil and Cervarix, composed of self-assembling L1 moieties in virus-like particle structures, have entered into market in 2006 and 2009, respectively (5). The Gardasil harbors a mix of HPV-16/18/6/11 L1-VLP, while Cervarix composed HPV-16/18 VLPs. Despite ideal results that were imagined for these approved vaccines, some drawback looks inevitable. One considerable shortcoming is their limited serotype coverage (9). The induced humoral immunity of these vaccines is primarily type specific, so this level of protection will not support complete protection from other oncogenic types. In sum, with respect to vaccine cost, restricted type coverage, and person-dependent efficacy of Gardasil and Cervarix for high-risk serotypes, designing a new preventive/therapeutic vaccine with extended coverage potency is an essential task (2,14).

The epitope-based vaccines can be designed for humoral and cellular stimulation and is very safe approach when compared to other strategies. B lymphocytes play a crucial role in the humoral immunity of adaptive immune system. The distinguishing of B cell epitopes as immunogenic regions within a target protein is the goal of some immunologic assays. Similar to cellular epitopes, B cell epitopes can play a vital role in the development of peptide vaccine (24). Identification of B cell epitopes reactive to HPV L1-derived peptides performed in animal model several times with the aims of designing effective epitope-based vaccine (7). The experimental techniques are time-consuming and labor-intensive, therefore several in silico methodologies are being created to distinguish epitopes of protein (6). On another side, computational techniques offer fast, easy, reliable, and cost-effective approaches for predicting B cell epitopes. Using bioinformatics tools, scientists can analyze the protein of interest for extraction of epitopes rather than potential binding sites in epitope-based vaccines (6,7). Furthermore, improved reliability of computational models in the prediction of desired B cell epitope will surely facilitate a preexperimental step of vaccine development.

Development of a new-generation peptide-based vaccine covering more carcinogenic serotypes may get rid of limited coverage of current HPV vaccines. Others also designed VLP-based vaccine by engineering of recombinant L1 VLP protein containing major preS1 epitopes of hepatitis B virus using bioinformatics analysis (10). Also, in silico analysis was applied to the development of peptide vaccine targeting different HPV proteins (11). They surveyed along 20 L1 sequence, and afterward, three conserved B cell peptides were extracted for next experimental assay. Recently, in an attempt to design multiepitope therapeutic vaccine against cancer, one group employed the benefit of computational models (15).

In this study, a multiepitope vaccine against major HPV structural proteins L1/L2 was designed by immunoinformatics models considering the coverage potency for HPV 16, 18, 31, 45, 6, and 11 subtypes as a universal vaccine. These proteins are considered to be immunogenic for vaccines design. For this purpose, to obtain surface and conserve B cell epitopes in silico by a convenient computational method, we extended informatics strategies on different virus genes and tried to identify consensus B cell epitopes suitable for prophylactic regimen. To our knowledge, reports of B cell epitope prediction for L1 and L2 of HPVs is limited; so we predicted the most probable immunogenic regions of L1 and L2 among high-risk HPV 16, 18, 31, 6, and 45 rather than type 11 serotypes based on BCPREDS defaulted model, and solvent accessibility of structure was extrapolated. For determination of solvent accessibility, the precise structure of protein is necessary. The crystal structure of L1 and L2 proteins was not determined. So we constructed a three-dimensional (3D) molecular model of the L1 by Swiss Model server, but this server was unable to provide model of L2 protein. Finally, based on best estimated epitope regions extracted from molecular modeling analysis, N-glycosylated sites, codon optimization, and other bioinformatics analyses, B cell epitope relevant to L1 and L2 protein was extracted and modified and then presented to new HPV constructs for the next experimental assessment.

Materials and Methods

Protein sequences collection

Protein reference sequences of L1 and L2 from common/high-risk HPVs, including type 6, 11, 16, 18, 31, and 45, are retrieved from http://ncbi.nlm.nih.gov/Databases. The sequence of L1 and L2 from HPV 6 and 11 was completely similar; thus, we omitted the sequences of L1 and L2 from HPV 6 for the next analysis and sequence of L1 and L2 from HPV 11 was considered for both 6 and 11 serotypes for further analysis. The sequences and corresponding accession number are listed in Table 1.

Table 1.

Sequences Used in This Study for Designing Construct

Serotype	L1 (Accession)	L2 (Accession)	Importance
HPV type 11	CCB84764	CCB84763	High prevalent/moderate risk
HPV type 16	AAV91691	AAO85414	High risk
HPV type 18	AAP20601	AGU90429	High risk
HPV type 31	AEI61093	AEI61092	High risk
HPV type 45	AGU90648	AGU90647	High risk

HPV, human papillomavirus.

B cell epitope prediction

For the prediction of epitope regions in L1 and L2 proteins, the BCPREDS server (http://ailab.cs.iastate.edu/bcpreds/predict.html) was employed based on recommended defaults. This server allows users to choose a modeling assay among three available prediction methods: (1) AAP method, (2) BCPred, and (3) FBCPred. In our prediction assay, BCPred method was selected because it could be possible to optionally specify the desired epitope length and specificity, its threshold following antigen sequence submission. Fixed length epitope was 20 mers peptide, specificity was 75% and used overlap filter in above servers.

Solvent accessibility of epitopes

After selection of epitopes, these epitopes were refined by the degree of solvent accessibility. Each residue at the surface of a protein can potentially be touched by water, and the area of an atom on the surface that can be touched by water is called the accessible molecular surface or solvent-exposed area. So it is advisable to select epitopes with higher degree of solvent accessibility. For fulfillment it, the structure of the L1 and L2 proteins was needs. The crystal structures of L1 and L2 protein had not been determined before. For obtaining a 3D structure model of L1 proteins, the crystal structure of these proteins was used as PDB template in Swiss Model Alignment interface protein modeling server (10). For L1 HPV 11, the crystal structure of L1 HPV 11 (PDB code: 2r5kD, Uniprot Entries Accession P04012); for L1 HPV 16, the crystal structure of L1 HPV 16 (PDB code: 1dzlA, Uniprot Entries Accession P03101); for L1 HPV 18, the crystal structure of L1 HPV 18 (PDB code: 2r5iG, Uniprot Entries Accession P06794); for L1 HPV 31, the crystal structure of L1 HPV 16 (PDB code: 1dzlA, Uniprot Entries Accession P03101); and for L1 HPV 45, the crystal structure of L1 HPV 18 (PDB code: 2r5iG, Uniprot Entries Accession P06794) was used. Despite submission sequence of L2 proteins, 3D structure models of L2 were not obtained. The fit between the 3D structures of models was evaluated in SWISS-PDB Viewer by calculating the root mean square deviation after iterative fitting. After preparation of structures, the solvent accessibility of these epitopes was studied in SPDBV software. Finally, the epitopes with a high degree of solvent accessibility were selected.

N-glycosylated sites

In N-glycosylation, carbohydrate chains are attached to nitrogen of asparagine and serve various functions. N-linked glycosylation is the most common type of glycosidic bond and is important for the folding of some eukaryotic proteins and occurs in eukaryotes, but very rarely in bacteria. We designed a final expression of construct in Escherichia coli, therefore it is necessary to delete the epitopes that have N-glycosylated sites. For prediction of N-glycosylated sites, the ensemblegly-Iowa server was used. Epitopes which selected from the previous step were analyzed in the ensemblegly-Iowa server and sequences that have an N-glycosylated site were excluded in the design of final construct.

Reverse translation of first structure

After selection of best epitopes of L1 and L2 based on bioinformatics analysis, these epitopes fused to each other, and the first structure of this construct was created by the Lasergene seqbuilder server. The final construct back translated into nucleotide sequences using reverse translate tool of Sequence Manipulation Suite (www.bioinformatics.org/sms2/rev_trans.html). The Sequence Manipulation Suite is a collection of JavaScript programs for analyzing the short DNA and protein sequences. In this server, reverse translate accepts a protein sequence as input and applies a codon usage table to generate a DNA sequence representing the most likely nondegenerate coding sequence. A consensus sequence derived from all the possible codons for each amino acid is also returned.

Codon optimization, rare codon determination

Codon optimization is a technique that supports the efficiency of DNA expression vectors used in DNA vaccination to achieve optimum expression of a foreign gene in the host's cell body (18). Constructs were optimized using the OPTIMIZER tool in genome server (http://genomes.urv.cat/OPTIMIZER), which evaluates the optimal codon corresponding to each amino acid encoded by the DNA and replaces the codons. OPTIMIZER is an online PHP application that optimizes the codon usage of a DNA sequence to increase its expression level. Users can introduce their own preference tables to be used in the optimization process or use precomputed tables from >150 prokaryotic species under a strong translational selection. This server can be useful for predicting and optimizing the level of expression of a gene in heterologous gene expression with the most optimal codon. Evaluation of final construct for rare codon analysis to enhance gene expression level and protein solubility has a very important step. For this purpose, we used rare codon analysis tools deposited in GenScript's server (www.Genscript.com/cgi-bin/tools/rare_codon_analysis). This tool reads protein coding DNA sequence and calculates its organism-related properties, such as CAI (Codon Adaptation Index), guanine/cytosine (GC) content, and protein codons and prime frequency distribution. This information is very useful to predict that the final construct gene will express well in E. coli. The CAI of construct improved to reach >0.8.

Design of vaccine constructs

Protein analyses such as surface probability plot, alpha and beta aliphatic region, translation of antigenic sequence, antigenicity index, and prediction of binding affinity of epitopes were measured as last check point for validity and reliability of construct by Hopp and Woods method in Expasy tools (http://web.expasy.org/protscale). At final step, pET28a vector (Invitrogen) selected as an expression plasmid for construction of multiepitope-based vaccine. The pET expression system is one of the most widely used systems for the cloning and in vitro expression of recombinant proteins in E. coli. In the final construct, the selected peptide sequence was inserted in BamHI/HindIII restriction sites. Furthermore the ATG initiation codon and GGATCC sequence (BamHI restriction site) were introduced to aminoterminal and AAGCTT sequence (HindIII restriction site) was added to carboxyl end of expressed protein in final construct. Only epitopes with the highest antigenicity score were selected, then they fused to each other, and the sequence was ordered. Then, the DNA fragments of construct, digested by BamHI/HindIII, were inserted into the BamHI/HindIII restriction sites of digested/dephosphorylated pET28a vector.

Results

B cell epitope prediction

By the use of BCPred method in BCPREDS server of epitope regions in L1 and L2, protein sequences were determined. A total of 71 epitopes for L1 and 70 epitopes for L2 HPV serotypes were predicted. In the next step, epitopes with high score were used as candidates for application in vaccine construct (Table 2).

Table 2.

Predicted Epitopes in HPV 11 and 16 Serotypes by Using BCPREDS Server

Serotype	Protein	Position	Epitope	Score
HPV 11	L1	128	NSGGYGGNPGQDNRVNVGMD	1
		4	PSDSTVYVPPPNPVSKVVAT	1
		482	RPAVSKPSTAPKRKRTKTKK	0.998
		300	EAQLFNKPYWLQKAQGHNNG	0.994
		169	TQCSNTSVQNGDCPPLELI	0.99
		422	ITCQKPTPEKEKQDPYKDMS	0.99
		400	FGLSPPPNGTLEDTYRYVQS	0.989
		99	CTGLEVGRGQPLGVGVSGHP	0.978
		275	KGGNNRSSVASSIYVHTPSG	0.972
		44	VGHPYYSIKKVNKTVVPKVS	0.959
		345	VSKSATYTNSDYKEYMRHVE	0.89
		227	KYPDYLQMAADPYGDRLFFY	0.887
		252	MFARHFFNRAGTVGEPVPDD	0.859
		190	SVIQDGDMVDTGFGAMNFAD	0.827
HPV 11	L2	71	IPLGSSPKPAITGGPAARPP	1
		197	FVVSSSDSGPTSSTPLPRAF	1
		105	SLIEESAIINAGAPEVVPPT	1
		159	FTEPSVIQPQPPVEASGHIL	0.999
		50	GVFFGGLGIGTGAGSGGRAG	0.999
		19	TCKATGTCPPDVIPKVEHTT	0.999
		403	PTASMGTPFSPVTPALPTGP	0.998
		363	VTQSYLTSTPNTLSQSWGNT	0.996
		240	STPQRLVTYDNPVYEGEDVS	0.974
		342	NDTFDIYAEPFDPIPDPVQH	0.957
		294	RFSRIGQRGSMYTRSGQHIG	0.948
		126	GGFTITSSESTTPAILDVSV	0.872
HPV 16	L1	511	RKATPTTSSTSTTAKRKKRK	1
		36	TVYLPPVPVSKVVSTDEYVA	1
		453	ACQKHTPPAPKEDPLKKYTF	1
		329	DAQIFNKPYWLQRAQGHNNG	0.999
		198	GSPCTNVAVNPGDCPPLELI	0.998
		426	EDWNFGLQPPPGGTLEDTYR	0.997
		72	GHPYFPIKKPNNNKILVPKV	0.994
		128	CVGVEVGRGQPLGVGISGHP	0.991
		152	LDDTENASAYAANAGVDNRE	0.983
		220	VIQDGDMVDTGFGAMDFTTL	0.964
		175	MDYKQTQLCLIGCKPPIGEH	0.921
		104	PDPNKFGFPDTSFYNPDTQR	0.919
		255	CKYPDYIKMVSEPYGDSLFF	0.856
		375	STSETTYKNTNFKEYLRHGE	0.803
HPV 16	L2	116	DAGAPTPVPSIPPDVSGFSI	1
		73	IPLGTRPPTATDTLAPVRPP	1
		204	STNPNTVTSSTPIPGSRPVA	1
		371	DDFITDTVTTPVPAIPSTSL	0.999
		344	ITPSTYTTTSHAASPTSINN	0.999
		416	IPINTTDQTPSLIPIVPGSP	0.999
		160	NPTFTDPSVLQPPTPAETGG	0.998
		137	TSTDTTPAILDINNTVTTVT	0.995
		21	TCKQAGTCPPDIIPKVEGKT	0.993
		51	MGVFFGGLGIGTGSGTGGRT	0.991
		237	VVDPAFVTAPTKLITYDNPA	0.975
		94	TVDPVGPSDPSIVSLVEETS	0.969
		393	YIPANTTIPFGGAYNIPLVS	0.964
		304	SRIGNKQTLRTRSGKSIGAK	0.961
		448	LHPSYYMLRKRRKRLPYFFS	0.868
		258	EGIDVDNTFYFPSNDNSINI	0.753
HPV 18	L1	547	RSAPSATTSSKPAKRVRVRA	1
		107	GNPYFRVPAGGGNKQDIPKV	1
		457	NSSILEDWNFGVPPPPTTSL	0.996
		365	SQLFNKPYWLHKAQGHNNGV	0.989
		75	PPPSVARVVNTDDYVTRTSI	0.987
		165	GVEIGRGQPLGVGLSGHPFY	0.98
		339	KGTGMRASPGSCVYSPSPSG	0.979
		407	ASTQSPVPGQYDATKFKQYS	0.976
		197	TSNVSEDVRDNVSVDYKQTQ	0.956
		318	ARHFWNRAGTMGDTVPQSLY	0.937
		257	EDGDMVDTGYGAMDFSTLQD	0.894
		235	ACKSRPLSQGDCPPLELKNT	0.88
		297	QMSADPYGDSMFFCLRREQL	0.867
		12	YHLLPLYGPLYHPRPLPLHS	0.836
		526	LGRKFLVQAGLRRKPTIGPR	0.818
		139	PDPNKFGLPDTSIYNPETQR	0.786
HPV 18	L2	82	VDVGPTRPPVVIEPVGPTDP	1
		120	PTFTGTSGFDITSAGTTTPA	1
		408	PSTTSVWPIVSPTAPASTQY	0.999
		20	TCKQSGTCPPDVVPKVEGTT	0.998
		351	YADDMDPAVPVPSRSTTSFA	0.997
		206	GEEPISSTPLPTVRRVAGPR	0.996
		58	GIGTGSGTGGRTGYIPLGGR	0.994
		387	VTVPLTSSWDVPVYTGPDIT	0.988
		153	STTNFTNPAFSDPSIIEVPQ	0.987
		319	FYHDISPIAPSPEYIELQPL	0.962
		438	WPLYYFIPKKRKRVPYFFAD	0.93
		248	LITYDNPAFEPVDTTLTFDP	0.888
		298	RLGQRATMFTRSGTQIGARV	0.813
HPV 31	L1	485	KRSAPSASTTTPAKRKKTKK	1
		13	LPPVPVSKVVSTDEYVTRTN	1
		128	DDTENSNRYAGGPGTDNREC	1
		172	KGSPCSNNAITPGDCPPLEL	0.995
		304	DAQIFNKPYWMQRAQGHNNG	0.995
		47	HPYYSIPKSDNPKKIVVPKV	0.992
		427	ITCQKTAPQKPKEDPFKDYV	0.992
		103	CVGLEVGRGQPLGVGISGHP	0.956
		193	NSVIQDGDMVDTGFGAMDF	0.955
		397	PAILEDWNFGLTTPPSGSLE	0.941
		149	SMDYKQTQLCLLGCKPPIGE	0.907
		79	PDPNKFGFPDTSFYNPETQR	0.893
		238	MVAEPYGDTLFFYLRREQMF	0.819
		336	DTTRSTNMSVCAAIANSDTT	0.785
		271	SVPTDLYIKGSGSTATLANS	0.74
HPV 31	L2	118	GAPAPIPHPPTTSGFDIATT	1
		200	TNNENITSSTPIPGVRRPAR	1
		84	EASIPIRPPVSIDPVGPLDP	1
		412	EHAPTQVFPFPLAPTTPQVS	0.997
		55	FGGLGIGSGSGTGGRTGYVP	0.995
		155	NPTFTDPSVLQPPTPAETSG	0.995
		332	SIEMQPLGASATTTSTLNDG	0.993
		362	TVDTPATHNVSPSTALQSTS	0.984
		391	TTVPLSTGFDIPIFSGPDVP	0.971
		21	TCKAAGTCPSDVIPKIEHTT	0.963
		311	ATIGARVHYYYDISSINPAG	0.954
		433	FVDGGDFYLHPSYYMLKRRR	0.842
		288	TSRRNTVRYSRLGNKQTLRT	0.814
		241	PKQLITYENPAYETVNAEES	0.808
HPV 45	L1	503	GLRRRPTIGPRKRPAASTST	1
		198	KGTLCKPAQLQPGDCPPLEL	0.998
		314	TPGSCVYSPSPSGSITTSDS	0.997
		425	NSSILENWNFGVPPPPTTSL	0.996
		375	ASTQNPVPNTYDPTKFKHYS	0.994
		127	WACVGMEIGRGQPLGIGLSG	0.984
		456	VTCQKDTTPPEKQDPYDKLK	0.984
		336	FNKPYWLHKAQGHNNGICWH	0.975
		165	VITQDVRDNVSVDYKQTQLC	0.934
		74	PYFRVVPSGAGNKQAVPKVS	0.929
		40	PPPSVARVVNTDDYVSRTSI	0.918
		219	NTIIEDGDMVDTGYGAMDFS	0.897
		263	QMSADPYGDSMFFCLRREQL	0.867
		288	WNRAGVMGDTVPTDLYIKG	0.826
HPV 45	L2	82	VDVGPTRPPVVIDPVGPTDP	1
		348	DVYADFPPPASTTPSTINKS	1
		55	GGLGIGTGSGSGGRTGYVPL	0.999
		412	TPMWPSTSPTNAATSTYIGI	0.997
		115	SGAPVPTFTGTSGFEITSSG	0.997
		203	SGSGTEPISSTPLPTVRRVA	0.997
		372	KYSLTMPSTAASSYSNVTVP	0.994
		153	SSTSFTNPAFSDPSIIEVPQ	0.979
		20	TCKQSGTCPPDVINKVEGTT	0.972
		433	GTQYYLWPWYYYFPKKRKRI	0.955
		304	TMFTRSGKQIGGRVHFYHDI	0.94
		249	VTFDNPAYEPLDTTLSFEPT	0.926

Solvent accessibility of epitopes

After position extraction, suggested epitopes were refined by the degree of solvent accessibility. To select epitopes with a higher degree of solvent accessibility, the structure of the L1 and L2 proteins was studied. The 3D structure models of L1 were obtained in the Swiss Model Alignment interface protein modeling server due to lack of the crystal structure for this protein. The results of homology modeling for L1 proteins from HPV types 11, 16, 18, 31, and 45 are listed in Table 3. As previously mentioned, this server was unable to model L2 protein. For further analysis of epitopes in L2 proteins, multiple sequence alignment of L2 protein from retrieved HPV 11, 16, 18, 31, and 45 serotype demonstrated that amino acid residue 11–200 have 50% homology. Therefore, extracted epitopes from this area were selected as epitopes of choice for the next step. Thereafter, the solvent accessibility of these epitopes was studied in the SPDBV software. Finally, as shown in Figure 1, the epitopes with a high degree of solvent accessibility were selected.

FIG. 1.

Solvent accessibility of epitopes. The residue main chain and side chains are represented in stick form (in white). Individual epitopes are identified and represented in Space-Filling model. The sequences of these epitopes are represented in Tables 2 and 5.

Table 3.

The Results of Homology Modeling and Model Information in SWISS-MODEL Workspace

	Modeled residue range	Template	Sequence identity [%]	Z-score
L1-HPV 11	20 to 470	2r5kD	92.24	−3.8
L1-HPV 16	46 to 500	1dzlA	99.12	−4
L1-HPV 18	82 to 536	2r5iG	92.09	−2.24
L1-HPV 11	20 to 474	1dzlA	83.96	−4.04
L1-HPV 11	47 to 504	2r5iG	81.44	−2.77

N-glycosylated sites

By the use of ensemblegly-Iowa server, N-glycosylation site of these epitopes was predicted. Epitopes analyzed in this server and position of asparagine that glycosylated were determined. To ensure for expression of construct in E. coli, epitopes that have an N-glycosylated site were excluded in the design of final construct. The score and potential of asparagine residue in the sequence of epitopes are demonstrated in Table 4.

Table 4.

N-Linked Glycosylated Sites

	Position	Residue	Score	Threshold	Prediction		Position	Residue	Score	Threshold	Prediction
A						B
	22	N	0.031	0.3	−		7	N	0.002	0.3	−
	29	N	0.012	0.3	−		142	N	0.015	0.3	−
	38	N	0.694	0.3	+		175	N	0.826	0.3	+
	60	N	0.002	0.3	−		184	N	0.001	0.3	−
	76	N	0.004	0.3	−		278	N	0.002	0.3	−
	100	N	0.808	0.3	+		297	N	0.011	0.3	−
	123	N	0.002	0.3	−		370	N	0.911	0.3	+
	166	N	0.003	0.3	−		401	N	0.041	0.3	−
	173	N	0.007	0.3	−		409	N	0.922	0.3	+
	180	N	0.001	0.3	−
	185	N	0.002	0.3	−
	188	N	0.006	0.3	−
	219	N	0.83	0.3	+
	224	N	0.044	0.3	−
	251	N	0.017	0.3	−
	258	N	0.834	0.3	+
	304	N	0.006	0.3	−
	323	N	0.01	0.3	−
	324	N	0.832	0.3	+
	350	N	0.004	0.3	−
	362	N	0.005	0.3	−
	363	N	0.014	0.3	−
	369	N	0.003	0.3	−
	383	N	0.993	0.3	+
	398	N	0.013	0.3	−
	436	N	0.795	0.3	+
	444	N	0.003	0.3	−
	452	N	0.77	0.3	+
	491	N	0.006	0.3	−
C						D
	17	N	0.015	0.3	−		7	N	0.002	0.3	−
	27	N	0.002	0.3	−		177	N	0.983	0.3	+
	49	N	0.01	0.3	−		178	N	0.03	0.3	−
	76	N	0.005	0.3	−		187	N	0.002	0.3	−
	79	N	0.007	0.3	−		188	N	0.989	0.3	+
	118	N	0.01	0.3	−		220	N	0.004	0.3	−
	142	N	0.003	0.3	−		234	N	0.003	0.3	−
	143	N	0.006	0.3	−		236	N	0.004	0.3	−
	144	N	0.001	0.3	−		282	N	0.001	0.3	−
	167	N	0.002	0.3	−		292	N	0.002	0.3	−
	178	N	0.003	0.3	−		299	N	0.011	0.3	−
	210	N	0.002	0.3	−		301	N	0.01	0.3	−
	217	N	0.886	0.3	+		304	N	0.001	0.3	−
	224	N	0.006	0.3	−		336	N	0.003	0.3	−
	229	N	0.003	0.3	−		362	N	0.002	0.3	−
	263	N	0.003	0.3	−		390	N	0.003	0.3	−
	267	N	0.002	0.3	−		391	N	0.008	0.3	−
	278	N	0.031	0.3	−		425	N	0.992	0.3	+
	302	N	0.937	0.3	+		435	N	0.004	0.3	−
	348	N	0.012	0.3	−		447	N	0.962	0.3	+
	356	N	0.005	0.3	−
	371	N	0.003	0.3	−
	376	N	0.006	0.3	−
	394	N	0.002	0.3	−
	406	N	0.003	0.3	−
	407	N	0.022	0.3	−
	413	N	0.003	0.3	−
	427	N	0.996	0.3	+
	443	N	0.011	0.3	−
	445	N	0.012	0.3	−
	481	N	0.971	0.3	+
	489	N	0.012	0.3	−
	536	N	0.003	0.3	−
E						F
	8	N	0.015	0.3	−		22	N	0.018	0.3	−
	14	N	0.002	0.3	−		28	N	0.01	0.3	−
	87	N	0.009	0.3	−		37	N	0.81	0.3	+
	89	N	0.006	0.3	−		123	N	0.007	0.3	−
	105	N	0.002	0.3	−		200	N	0.998	0.3	+
	119	N	0.041	0.3	−		203	N	0.004	0.3	−
	143	N	0.005	0.3	−		223	N	0.01	0.3	−
	154	N	0.003	0.3	−		282	N	0.002	0.3	−
	177	N	0.002	0.3	−		297	N	0.004	0.3	−
	188	N	0.002	0.3	−		389	N	0.001	0.3	−
	220	N	0.003	0.3	−		430	N	0.998	0.3	+
	234	N	0.994	0.3	+
	242	N	0.515	0.3	+
	288	N	0.028	0.3	−
	358	N	0.018	0.3	−
	404	N	0.004	0.3	−
	416	N	0.004	0.3	−
	417	N	0.002	0.3	−
	423	N	0.997	0.3	−
	437	N	0.937	0.3	+
	492	N	0.74	0.3	+
	500	N	0.003	0.3	−
	534	N	0.002	0.3	−
	545	N		0.3	−
G						H
	3	N	0.01	0.3	−		10	N	0.013	0.3	−
	30	N	0.018	0.3	−		18	N	0.003	0.3	−
	64	N	0.007	0.3	−		194	N	0.938	0.3	+
	89	N	0.003	0.3	−		226	N	0.005	0.3	−
	114	N	0.002	0.3	−		240	N	0.006	0.3	−
	125	N	0.002	0.3	−		241	N	0.003	0.3	−
	157	N	0.003	0.3	−		243	N	0.997	0.3	+
	164	N	0.011	0.3	−		288	N	0.003	0.3	−
	166	N	0.087	0.3	−		295	N	0.01	0.3	−
	176	N	0.003	0.3	−		304	N	0.892	0.3	+
	210	N	0.004	0.3	−		308	N	0.002	0.3	−
	211	N	0.011	0.3	−		331	N	0.009	0.3	−
	225	N	0.039	0.3	−		340	N	0.002	0.3	−
	252	N	0.004	0.3	−		366	N	0.001	0.3	−
	259	N	0.005	0.3	−		388	N	0.006	0.3	−
	295	N	0.882	0.3	+		409	N	0.809	0.3	+
	321	N	0.951	0.3	+		429	N	0.947	0.3	+
	341	N	0.003	0.3	−
	353	N	0.003	0.3	−
	354	N	0.02	0.3	−
	360	N	0.003	0.3	−
	374	N	0.991	0.3	+
	383	N	0.007	0.3	−
	392	N	0.004	0.3	−
	428	N	0.002	0.3	−
	436	N	0.01	0.3	−
	483	N	0.004	0.3	−
I						J
	7	N	0.002	0.3	−		7	N	0.002	0.3	−
	32	N	0.004	0.3	−		61	N	0.006	0.3	−
	45	N	0.007	0.3	−		107	N	0.007	0.3	−
	47	N	0.009	0.3	−		187	N	0.004	0.3	−
	77	N	0.032	0.3	−		207	N	0.008	0.3	−
	101	N	0.004	0.3	−		259	N	0.007	0.3	−
	113	N	0.002	0.3	−		281	N	0.002	0.3	−
	136	N	0.005	0.3	−		298	N	0.019	0.3	−
	147	N	0.001	0.3	−		393	N	0.752	0.3	+
	179	N	0.004	0.3	−		415	N	0.994	0.3	+
	201	N	0.637	0.3	+		450	N	0.008	0.3	−
	247	N	0.024	0.3	−
	317	N	0.016	0.3	−
	338	N	0.007	0.3	−
	365	N	0.002	0.3	−
	377	N	0.004	0.3	−
	378	N	0.015	0.3	−
	384	N	0.003	0.3	−
	398	N	0.998	0.3	+
	407	N	0.001	0.3	−
	411	N	0.006	0.3	−
	453	N	0.97	0.3	+
	459	N	0.01	0.3	−
	461	N	0.189	0.3	−

The residues that accept a glycosidic bond is shown with a plus sign (A: L1 HPV 11, B: L2 HPV 11, C: L1 HPV 16, D: L2 HPV 16, E: L1 HPV 18, F: L2 HPV 18, G: L1HPV 31, H: L2 HPV 31, I: L1 HPV 45, J: L2 HPV 45).

Selecting epitopes and fusion of peptides for generating vaccine constructs

After the evaluation of epitope scores, two optimal epitopes were selected for each HPV type. L1 and L2 protein epitopes of serotype 6 and 11 were similar; therefore, we selected one epitope for both. As shown in Table 5, finally 20 epitopes were selected. These epitopes were connected to each other and the first structure of this construct was created (Fig. 2). The tandem order of them was determined just by the order of their appearance in Table 5 and they are designated as E1 to E10 for L1 protein and E11 to E20 for L2 protein.

FIG. 2.

Design of primary structure for final construct from selected epitope. Number of HPV are determined in Table 5. HPV, human papillomavirus.

Table 5.

The Selected Epitopes from L1 and L2 Epitope Pool

No.	Serotype	Protein	Start	Sequence peptide
01	HPV 11	L1	128	NSGGYGGNPGQDNRVNVGMD
02	HPV 11	L1	169	GTQCSNTSVQNGDCPPLELI
03	HPV 11	L2	19	TCKATGTCPPDVIPKVEHTT
04	HPV 11	L2	105	SLIEESAIINAGAPEVVPPT
05	HPV 16	L1	152	LDDTENASAYAANAGVDNRE
06	HPV 16	L1	375	STSETTYKNTNFKEYLRHGE
07	HPV 16	L2	21	TCKQAGTCPPDIIPKVEGKT
08	HPV 16	L2	116	DAGAPTPVPSIPPDVSGFSI
09	HPV 18	L1	407	ASTQSPVPGQYDATKFKQYS
10	HPV 18	L1	457	NSSILEDWNFGVPPPPTTSL
11	HPV 18	L2	20	TCKQSGTCPPDVVPKVEGTT
12	HPV 18	L2	82	VDVGPTRPPVVIEPVGPTDP
13	HPV 31	L1	47	HPYYSIPKSDNPKKIVVPKV
14	HPV 31	L1	172	KGSPCSNNAITPGDCPPLEL
15	HPV 31	L2	21	TCKAAGTCPSDVIPKIEHTT
16	HPV 31	L2	84	EASIPIRPPVSIDPVGPLDP
17	HPV 45	L1	375	ASTQNPVPNTYDPTKFKHYS
18	HPV 45	L1	425	NSSILENWNFGVPPPPTTSL
19	HPV 45	L2	20	TCKQSGTCPPDVINKVEGTT
20	HPV 45	L2	82	VDVGPTRPPVVIDPVGPTDP

Evaluation of protein sequence

The sequences of primary construct were analyzed in Expasy's ProtParam tool and results are presented in Table 6. The results indicated that molecular weight of this protein was 41.8 kDa and has 400 amino acid residues. Because the total number of positively charged residues was nearly twice that of positively charged residues, the computed isoelectric point (pI) value was below 7 (4.71), which indicates the acidic nature of final product. The aliphatic index shows that this construct has low thermal stability. The Grand Average Hydropathicity (GRAVY) value was negative (−0.494). Besides, amino acid composition revealed that most abundant amino acid residues in our construct were glutamine and leucine, respectively. The antigenicity of construct was analyzed by protean program. This plot showed that major parts of new protein construct have hydrophilic property and thus harbor antigenic potency (Fig. 3).

FIG. 3.

The diagram of antigenicity determined by protean software for final construct.

Table 6.

Results of Designed Protein Characterization and Calculated Parameters from Expasy's ProtParam Tool

Parameters	E2 protein
No. of amino acids	400
Theoretical pI	4.71
Molecular weight	41.8597
Positively charged residues (Arg+Lys)	28
Negatively charged residues (Asp+Glu)	44
Formula	C1831H2875N485O606S15
Instability index	41.77
Total no. of atoms	5,812
AI, GRAVY	64.00, −0.494

AI, aliphatic index; GRAVY, grand average hydropathicity; pI, theoretical isoelectric point.

Reverse translation of first structure

After confirmation of the new protein antigenicity, amino acid sequences of this construct back translated into nucleotide by sequence manipulation suite server. This process was conducted based on codon usage table, and consensus sequence derived from all the possible codons is also returned. Reverse-translated sequence was evaluated for codon optimization and rare codon to enhance gene expression level and protein solubility. For this purpose, OPTIMIZER and GenScript's servers were employed together. This information is very useful to predict that the final construct gene will express well in E. coli host. A CAI of 1.0 is considered ideal, while a CAI of >0.8 is rated as good for expression in the desired expression organism. The lower the number contain the higher chance to gene be expressed poorly. GenScript's Optimum Gene™ codon optimization tool can typically improve the sequence to reach a CAI of >0.8. The GC content of gene was 54.26%. The ideal percentage range of GC content is between 30% and 70%. Any peaks outside this range will adversely affect transcriptional and translational efficiency. The percentage of low-frequency codons (<30%) based on our target host organism is 3%. This unoptimized gene employs tandem rare codons that can reduce the efficiency of translation or even disengage the translational machinery. Results are shown in Figures 4 and 5.

FIG. 4.

Codon optimization of final construct. Plot (A) shows relative addictiveness (wij) of construct before optimization (adaptation) to Escherichia coli codons, and plot (B) shows relative addictiveness (wij) of construct after optimization. The gray line shows the codon usage for each codon, black line depicts the mean codon usage, which was calculated from the E. coli. CAI value of the sequence before and after optimization was 0.63 and 0.97 and the GC content before and after optimization was 60.5 and 55.9, respectively. CAI, Codon Adaptation Index; GC, guanine/cytosine.

FIG. 5.

Representation of reverse translation and codon optimization. This screenshot shows translation of construct and optimization (adaptation) to E. coli codons.

Discussion

So far, three prophylactic HPV vaccines, based on VLP technology, have been approved and licensed. They are highly effective and immunogenic if administered before HPV exposure, which is to say, before the commencement of an active sexual life (18). L1-based prophylactic vaccines, such as Gardasil^®, Cervarix^®, and Gardasil 9^®, have been already licensed, while L2-based second-generation preventive vaccines are still under clinical trials. An important limitation of L1-based vaccines is that they induce type-restricted immune responses. L2-based vaccines can overcome this issue (18). The replication process of some viruses is error prone, so evolving of different strains is unavoidable; therefore, developing an effective vaccine with all serotypes coverage is nearly impossible. Although establishment of a universal vaccine for high evolving viruses such as influenza and HIV has many difficulties, some achievements have been published recently (16,23). Even providing an effective vaccine for HPV did not decrease the rate of cervical cancer in developing countries. HPV prevention remained a big concern, considering some shortages such as the presence of different genotypes and vaccine cost and availability (2,21). While around 15 high-risk HPVs were identified, current vaccines just protect humans from three types. The narrow coverage of current vaccine seems to be the most important limitation (2). The Gardasil and Cervarix harbor a mix of HPV-16/18/6/11 and HPV-16/18 L1 VLPs, respectively. Among all candidates, protein-based vaccine is a safe method for prevention and therapy (22). Development of low-cost and universal vaccine covering more carcinogenic serotypes may get rid of current HPV vaccines shortcoming (26). Theoretically, development of multiepitope-based vaccine containing common exposed sequences of HPV proteins induced wide spectrum of B cell immunity that consequently will establish potent humoral immunity. In fact, B lymphocytes play a crucial role in the humoral immunity of adaptive immune system. Similar to HLA class I and II epitopes, B cell epitopes also play a vital role in the development of protein vaccines. The B cell detects exposed parts of antigen and triggers potent humoral immune responses. The identification of exposed immunogenic regions within a protein, specifically those that elicit humoral immune response had been the goal of many studies before (2,25). Identification of B cell epitopes for HPV proteins was performed in animal models several times to design multiepitope-based vaccine (7,23). Identification of exposed motifs from virus antigens that contributed in binding to antibodies is helpful for this purpose (7,10). Immunobioinformatics analysis was applied for this kind of prediction. While experimental distinguishing of B cell epitopes is costly and time-consuming, in silico screening provided a convenient surrogate. In a previous study, Subramanian and Chinnappan (24) employed immunoinformatics tools to extract linear B cell epitope for E6 protein of high-risk HPVs. Kaushik et al. (12) by in silico analysis constructed peptide-based vaccine targeting L1 protein. In this study, we purposed a single polytopic candidate construct with six significant serotypes coverage (four high-risk HPVs and two high-rate HPV serotypes) as a new and universal vaccine candidate. Four high-risk and two low-risk HPV L1 and L2 protein sequences were evaluated for epitope predictions by the help of different immunoinformatics methods. At the end, 20 epitopes were selected as candidates for designing final construct and these peptides fused together in tandem without linker. After first epitope screening, nearly 140 epitopes were selected from 6 common HPV types. For the next step, 20 epitopes with higher degree of solvent accessibility and also higher degree of homology were selected. One of the most cited features of the epitope is hydrophilic property of surrounding sequence. This feature was first described by Novotny who calculated the solvent accessible surface area of epitopic residues (5). To select epitopes with higher degree of solvent accessibility, the 3D structure of the L1 and L2 proteins was constructed. For L2 protein, multiple alignments were performed on HPV 11, 16, 18, 31, and 45 serotypes to distinguish homology. It demonstrated that 50% homology among amino acid residue 11–200 present therefore extracted epitopes from this area selected as epitopes of choice for next step. These epitope sequences fused together in a tandem. Evaluation of new recombinant protein sequence indicated a molecular weight of 41.8 kDa with 400 amino acids. The computed isoelectric point (pI) value was 4.71, which indicates the acidic nature of final product. The aliphatic index showed that this construct has low thermal stability. Analyses of the GRAVY values showed a negative score of −0.494, indicative of suitable antigen solubility. Check of antigenicity also revealed that major parts of new protein construct have hydrophilic property and thus harbor antigenic potency. This finding also confirmed the exposure potency of almost all parts of newly constructed fusion polypeptide. Afterward, sequence of final construct was reverse translated to DNA sequence and then codon optimized. This codon-optimized sequence showed CAI of >0.8 for expression in E. coli and this score is acceptable. Finally, this sequence ligated into a pET28a bacterial expression vector through BamHI and HindIII restriction sites. Based on rational design of this construct, we expected that this construct would be a suitable candidate for universal vaccines against six common HPV types that are responsible for >90% of cervical cancer. Future generation of new highly immunogenic and effective vaccines such as L1/L2-based, capsomer-based, and pseudovirion-based vaccines can be further developed, thanks to computer-aided design and bioinformatics/computational biology (18). This kind of prediction and design by the help of bioinformatics will alleviate cost, time, labor, and space of animal experiments, especially when constructing a vaccine. Our designed construct unenviably need further experimental confirmation but at least validated as a universal derived antigen when evaluated computationally. Next, the universal vaccine will be expressed and immune responses against vaccine candidate would be evaluated in an experimental animal model.

Footnotes

Author Disclosure Statement

No competing financial interests exist. None of the work was paid for by the Iranian government or performed by its employees.

References

Arbyn

, Castellsague

, de Sanjose

, et al. Worldwide burden of cervical cancer in 2008. Ann Oncol, 2011; 22:2675–2686.

Barrios

, and Celis

. TriVax-HPV: an improved peptide-based therapeutic vaccination strategy against human papillomavirus-induced cancers. Cancer Immunol Immunother, 2012; 61:1307–1317.

Bharadwaj

, Hussain

, Nasare

, and Das

. HPV & HPV vaccination: issues in developing countries. Indian J Med Res, 2009; 130:327–333.

Doorbar

, Quint

, Banks

, et al. The biology and life-cycle of human papillomaviruses. Vaccine, 2012; 30 Suppl 5:F55–F70.

Einstein

, Baron

, Levin

, et al. Comparison of the immunogenicity and safety of Cervarix and Gardasil human papillomavirus (HPV) cervical cancer vaccines in healthy women aged 18–45 years. Human Vaccin, 2009; 5:705–719.

Eisenberg

, Wilcox

, and McLachlan

. Hydrophobicity and amphiphilicity in protein structure. J Cell Biochem, 1986; 31:11–17.

El-Manzalawy

, and Honavar

. Recent advances in B-cell epitope prediction methods. Immunome Res, 2010; 6 Suppl 2:S2.

Fukui

, Matsueda

, Kawano

, et al. Identification of B cell epitopes reactive to human papillomavirus type-16L1-derived peptides. Virology J, 2012; 9:199.

Guan

, Howell-Jones

, Li

, et al. Human papillomavirus types in 115,789 HPV-positive women: a meta-analysis from cervical infection to cancer. Int J Cancer, 2012; 131:2349–2359.

10.

Jagu

, Karanam

, Gambhira

, et al. Concatenated multitype L2 fusion proteins as candidate prophylactic pan-human papillomavirus vaccines. J Natl Cancer Inst, 2009; 101:782–792.

11.

Joshi

, Lewis

, Singharoy

, et al. Epitope engineering and molecular metrics of immunogenicity: a computational approach to VLP-based vaccine design. Vaccine, 2013; 31:4841–4847.

12.

Kaushik

, Kumar Shrivastav

, Shrivastav

, et al. In silico vaccine design against the target L1 binding protein of human papillomavirus, an etiological agent of cervical cancer using bioinformatics tools. Int J Pharm Sci Res, 2013; 4:4758–4763.

13.

Knipe

, and Howley

. 2013. Fields Virology. 6th, ed. Wolters Kluwer, Philadelphia.

14.

, Maraj

, Tran

, et al. Emerging human papillomavirus vaccines. Expert Opin Emerg Drugs, 2012; 17:469–492.

15.

, Roden

, and Wu

. Current status of human papillomavirus vaccines. J Formos Med Assoc, 2010; 109:481–483.

16.

Nezafat

, Ghasemi

, Rahbar

, et al. A novel multi-epitope peptide vaccine against cancer: an in silico approach. J Theor Biol, 2014; 21:121–134.

17.

Novotny

, Handschumacher

, Haber

, et al. Antigenic determinant in proteins coincide with surface regions accessible to large probes (antibody domains). Proc Nat Acad Sci U S A, 1986; 83:226–230.

18.

Panatto

, Amicizia

, Bragazzi

, et al. Human papillomavirus vaccine: state of the art and future perspectives. Adv Protein Chem Struct Biol, 2015; 101:231–322.

19.

Pica

, and Palese

. Toward a universal influenza virus vaccine: prospects and challenges. Ann Rev Med, 2013; 64:189–202.

20.

Puigbo

, Guzman

, Romeu

, et al. OPTIMIZER: a web server for optimizing the codon usage of DNA sequences. Nucleic Acids Res, 2007; 35(Web Server issue):W126–W131.

21.

Schiffman

, Castle

, Jeronimo

, et al. Human papillomavirus and cervical cancer. Lancet, 2007; 370:890–907.

22.

Schwede

, Kopp

, Guex

, et al. SWISS-MODEL: an automated protein homology-modeling server. Nucleic Acids Res, 2003; 31:3381–3385.

23.

Selvey

, Tindle

, Geysen

, et al. Identification of B-epitopes in the human papillomavirus 18 E7 open reading frame protein. J Immunol, 1990; 145:3105–3110.

24.

Subramanian

, and Chinnappan

. Prediction of promiscuous epitopes in the e6 protein of three high risk human papilloma viruses: a computational approach. Asian Pac J Cancer Prev, 2013; 14:4167–4175.

25.

Sun

, Ju

, Liu

, et al. Bioinformatics resources and tools for conformational B-cell epitope prediction. Comput Math Methods Med, 2013; 2013:943636.

26.

Tumban

, Peabody

, et al. A universal virus-like particle-based vaccine for human papillomavirus: longevity of protection and role of endogenous and exogenous adjuvants. Vaccine, 2013; 31:4647–4654.