Genetic variability in minor capsid protein (L2 gene) of human papillomavirus type 16 among Indian women.

Human papillomavirus type 16 (HPV-16) is the predominant genotype worldwide associated with invasive cervical cancer and hence remains as the focus for diagnostic development and vaccine research. L2, the minor capsid protein forms the packaging unit for the HPV genome along with the L1 protein and is primarily associated with transport of genomic DNA to the nucleus. Unlike L1, L2 is known to elicit cross-neutralizing antibodies and thus becomes a suitable candidate for pan-HPV prophylactic vaccine development. In the present study, a total of 148 cervical HPV-16 isolates from Indian women were analyzed by PCR-directed sequencing, phylogenetic analysis and in silico immunoinformatics tools to determine the L2 variations that may impact the immune response and oncogenesis. Ninety-one SNPs translating to 35 non-synonymous amino acid substitutions were observed, of these 16 substitutions are reported in the Indian isolates for the first time. T245A, L266F, S378V and S384A substitutions were significantly associated with high-grade cervical neoplastic status. Multiple substitutions were observed in samples from high-grade cervical neoplastic status as compared to those from normal cervical status (p = 0.027), specifically from the D3 sub-lineage. It was observed that substitution T85A was part of both, B and T cell epitopes recognized by MHC-I molecules; T245A was common to B and T cell epitopes recognized by MHC-II molecules and S122P/A was common to the region recognized by both MHC-I and MHC-II molecules. These findings reporting L2 protein substitutions have implications on cervical oncogenesis and design of next-generation L2-based HPV vaccines.

Introduction

Cervical cancer is the fourth common cancer occurring in women worldwide, with an estimated 570,000 new cases and 250,000 deaths occurring annually [1, 2]. Infection with high-risk or oncogenic human papillomavirus (HPV) types plays a critical role in the development of cervical cancer. HPV type 16 (HPV-16) remains the most prevalent high-risk type globally and in India [2, 3]. The ~8 kb genome of HPV-16 is packaged in the viral capsid which is composed of 360 molecules of the major capsid protein L1 and up to 72 molecules of the minor capsid protein L2, which largely lies buried inside L1 [4, 5]. The capsid proteins can assemble into virus-like particles (VLPs) and induce highly neutralizing antibodies that have shown promise in the prevention of papillomavirus-associated cancers. The currently available HPV vaccines, Gardasil and Cervarix are based on the L1, major capsid protein. However, targeting the L1 protein alone, prevents infection by limited HPV types because of lack of cross-reactive epitopes toward other HPV types. The HPV L2 protein on the other hand offers a distinct advantage, as it not only elicits neutralizing antibodies but also cross-neutralizing antibodies for different HPV types, especially against its ‘N’ terminus [4, 6, 7]. Therefore, development of a pan-HPV prophylactic vaccine derived from L2 regions seems feasible and needs to be explored.

For the development of an efficient L2 vaccine construct, recognition of appropriate epitope sequence is imperative to generate an efficient immune response [8, 9]. Implementation of informatics tools and online databases can facilitate this and reduce the time and experimental costs involved [10]. Determination of genetic variability, pertaining to epitopes, is important, as substitution of one or more amino acid(s) in the L2 protein could impact the conformation of epitopes relevant for viral neutralization. L2 protein also plays an important role during HPV infection, by binding to a secondary viral receptor to facilitate exit from the endosomes and delivery of the viral genome to appropriate domains within the nucleus [11]. Variations in L2 gene might lead to amino acid substitutions and the associated conformational changes can impact viral assembly and clinical outcome of infection. Ideally, vaccine constructs and diagnostics need to be developed locally and hence it is important to understand the geographical variations.

Data with regard to molecular analysis of L2 protein of HPV-16 from India is negligible. The present study was thus carried out to determine the genetic variability in the L2 protein of HPV-16 and to analyze the association of the identified variations with cervical disease status. We further predicted the immunogenic B and T cell epitopes of L2 using in silico immunoinformatics tools and identified the substitutions that may impact these epitopes.

Methods

The study was approved by Institutional Ethics Committee of ICMR-National AIDS Research Institute (ICMR-NARI), Pune, India [NARI/EC/approval/17-18/31].

Clinical specimens

A total of 148 cervical specimens positive for HPV-16 infection collected from women undergoing cervical screening from different geographical locales from India and stored at −70 °C as part of the sample repository in the Microbiology laboratory, ICMR-NARI were included in the study. The samples were classified as belonging to normal (n = 41), low-grade (n = 45) or high-grade (n = 62) cervical status as described previously [12]. HPV detection and genotyping for the samples was done using the Linear Array HPV genotyping test (Roche Molecular Systems, USA).

Nucleic acid extraction and sequencing

DNA was extracted from cervical samples using the QIAamp DNA mini kit (Qiagen, USA) and analyzed for L2 (3373–4794 bp) using two pairs of gene specific primers [HPV-16 L2 F1- TTACTTAACAATGCGACACA, HPV-16 L2 R1-TTATCCACATCTATACCTTCA, HPV-16 L2 F2-CCCTGCTTTTGTAACCACTC and HPV-16 R2-CGTGCAACATATTCATCCGT [13]. The amplified PCR products were purified with QIAquick PCR purification kit (Qiagen, USA) and sequenced using Big Dye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, USA). The sequence data was curated using SeqScape (V.2.6) and aligned in SEED 2 software [14]. MEGA7 [15] was used to identify site specific nucleotide variations and their corresponding amino acid substitutions using the reference entry, Refseq_ID: NC_001526.4 [16].

Phylogenetic analysis

Phylogenetic analysis of the 148 isolates of L2 gene sequences generated in this study was carried out using both alignment-based Maximum Likelihood (ML) method (for 1000 bootstrap runs) implemented in IQTREE server (http://iqtree.cibiv.univie.ac.at/) [17] and Return Time Distribution (RTD) (http://bioinfo.unipune.ac.in/RTD/home.htm) developed in-house which is an alignment-free method [18] to analyze their clustering pattern. A dataset consisting of total 211 entries of L2 was used which includes 63 entries from GenBank with known lineages (reference dataset; https://www.ncbi.nlm.nih.gov/nucleotide/) and 148 sequences generated in this study. Typing of the HPV-16 Indian isolates was done on basis of their L1 gene [12] and these lineage and sub-lineages were assigned to the L2 sequences of corresponding isolates. The OTU (operational taxonomic unit) labels in the tree were generated as ‘lineage/sub-lineage_isolate-ID_gene’. The OTU labels of isolates for which lineage/sub-lineage could not be assigned due to non-availability of their corresponding L1 sequences were generated as ‘XX_ isolate-ID_gene’. Phylogenetic trees were visualized using the iTOL server (https://itol.embl.de/) [19]. Recombination detection was carried out using the RDP4 software [20], while selection pressure analysis was carried out using Datamonkey server (http://www.datamonkey.org/) and a stringent cut-off of 10–5 was applied.

B cell and T cell epitope prediction

Experimentally characterized B and T cell epitopes of L2, were extracted from Immune Epitope Database (IEDB, http://www.iedb.org/home_v3.php). We further predicted linear B cell epitopes for L2 protein using Bepipred Linear Epitope Prediction 2.0 [21], Chou & Fasman Beta-Turn Prediction [22], Emini Surface Accessibility Prediction [23], Karplus & Schulz Flexibility Prediction [24], Kolaskar & Tongaonkar Antigenicity [25] and Parker Hydrophilicity Prediction [26] methods available on IEDB resource. Consensus peptides predicted using all the six prediction methods were taken and overlapping antigenic regions were concatenated. T cell epitopes having potential to recognize both, MHC-I and MHC-II were predicted. MHC-I binding epitopes were predicted using MHC-I binding consensus [27] and MHC-I processing (Proteosome, TAP) methods [28]. The predicted epitopes were filtered using percentile rank, predicted IC50 < 50 nM, total score (combined score of proteasome, TAP and cleavage scores) and checked for consensus. MHC-II binding epitopes were predicted using MHC-II binding consensus method [29, 30] and filtered using both, percentile rank and predicted IC50 (<50 nM). The overlapping MHC-I and MHC-II peptides were concatenated to yield non-redundant antigenic regions.

Statistical analysis

Statistical analysis was done using SPSS (V.15.0). HPV-16 L2 variations and cervical lesion grade were compared using Mann–Whitney U test as described earlier [12]. Odds ratio and 95% confidence intervals were computed to determine the association between normal and high-grade cervical disease status. Differences were considered to be statistically significant if p < 0.05.

Results

HPV-16 L2 variations

A total of 148 HPV-16 positive cervical samples were analyzed to study the extent of variations in L2 gene that resulted in detection of 91 nucleotide polymorphisms, of which, 43 (47.3%) were non-synonymous and 48 (52.7%) were synonymous nucleotide variations, translating to 35 (39.8%) non-synonymous and 53 (60.2%) synonymous amino acid substitutions, with no indels.

A total of 16 amino acid substitutions L75F, T85A, T94A, S122A, S134R, S270N, D272N, N273S, D334T, Q342L, T352P, H354Q, T377S, L390F, I418M, S426A in L2 protein are reported in the Indian HPV-16 isolates for the first time. Amino acid substitutions L330F (75.6%), S269P (28.4%) and D334N (24.3%) were most frequently observed.

Evolutionary analysis of HPV-16 L2

Phylogenetic trees depict the evolution of L2 gene of Indian isolates (Fig. 1). Similar clustering patterns were observed for both, alignment-based ML and alignment-free RTD methods (Fig. 1a, b) barring relative order of isolates as part of lineage assignments. Of the 148 Indian isolates, 134 isolates were assigned lineage on the basis of L1 gene, whereas 14 isolates could not be assigned lineage due to lack of corresponding L1 sequences. As can be seen in Fig. 1, 124/134 (92.5%) isolates clustered in accordance with lineage assignment based on L1. Thus, 107/124 (86.3%) and 17/124 (13.7%) isolates clustered with members of lineage A and lineage D, respectively. Of the 10/134 (7.5%) isolates, that did not cluster in accordance with their respective assigned lineages, 7 and 2 isolates with assigned lineage A clustered with members of lineage D and C respectively. The remaining 1 isolate with assigned lineage D clustered with members of lineage C. The 14 isolates without assigned lineage, clustered into A and D lineage clusters with 7 isolates in each. No evidence of recombination and positive selection was observed.

Fig. 1

Phylogenetic tree of HPV-16 L2 gene: a Alignment-based Maximum Likelihood tree and b Alignment-free RTD-based tree

HPV-16 L2 variations and cervical disease status

The association of amino acid substitutions with cervical lesion grade is shown in Table 1. Significantly higher occurrence of multiple substitutions in samples with high-grade cervical status as compared to normal cervical status (p = 0.027) was observed. T245A, L266F, S378V and S384A substitutions showed significant association with high-grade cervical status as compared to normal cervical status. These substitutions were predominantly observed in the D3 sub-lineage. The T245A substitution was observed in 19 samples (A1:2, A2:4, D1:1, D3:6 and unassigned sub-lineage:6), L266F was observed in 20 samples (A1:2, A2:4, D1:1, D3:7 and unassigned sub-lineage:6), S378V was observed in 20 samples (A1:2, A2:4, D1:1, D3:7 and unassigned sub-lineage:6) and S384A was observed in 20 samples (A1:2, A2:4, D1:1, D3:7 and unassigned sub-lineage:6). These 4 substitutions were found to be co-mutating in 19 samples.

Table 1

Human papillomavirus 16 amino acid substitutions in L2 by cervical lesion grade

HPV-16 L2 Amino acid substitutions	Cervical disease status			Total (n = 148)	ORa(95% CI)	p-value
	Normal (n = 41)	Low-grade (n = 45)	High-grade (n = 62)
L75F	2	0	0	2	0.73 (0.29, 1.78)	0.512
T85A	2	0	0	2	–	–
T94A	0	2	0	2	–	–
S122P	8	5	20	36	0.51 (0.19, 1.30)	0.180
S122A	0	3	0	3	–	–
S134R	1	0	0	1	–	–
T245A	2	3	19	24	0.12 (0.02, 0.53)	<0.01
L266F	2	4	20	26	0.10 (0.02, 0.49)	<0.01
S269P	10	9	23	42	0.55 (0.23, 1.32)	0.201
S270N	0	3	0	3	–	–
D272N	0	3	0	3	–	–
N273S	0	2	1	3	–	–
I306L	0	2	1	3	–	–
L330F	30	40	42	112	1.96 (0.77, 5.01)	0.180
T332S	1	0	9	10	0.15 (0.02, 1.21)	0.084
D334N	8	10	19	37	0.47 (0.18, 1.21)	0.124
D334H	0	0	1	1	–	–
D334T	0	0	1	1	–	–
E338D	0	5	3	8	–	–
Q342L	2	0	0	2	–	–
T352P	3	0	0	3	–	–
H354Q	0	0	2	2	–	–
T377S	0	0	1	1	–	–
S378V	5	5	20	30	0.29 (0.09, 0.85)	0.033
S378F	8	8	20	36	0.51 (0.19, 1.30)	0.180
S384A	5	5	20	30	0.29 (0.09, 0.85)	0.033
V385I	8	5	22	35	0.44 (0.17, 1.11)	0.120
L390F	0	0	1	1	–	–
I418M	0	0	1	1	–	–
I420T	8	5	20	33	0.51 (0.19, 1.30)	0.180
Q423H	0	3	1	4	–	–
A424T	8	5	20	33	0.51 (0.19, 1.30)	0.180
S426A	2	0	0	2	–	–
I428L	2	5	5	12	0.58 (0.11, 3.16)	0.699
A443G	7	5	20	32	0.43 (0.16, 1.14)	0.110

OR Odds ratio, 95% CI 95% confidence interval

aOR with 95% CI are presented for normal versus high-grade cervical status. The p values of significant associations are marked in boldface

Epitope prediction

Six experimentally validated non-overlapping B cell epitopes were obtained from IEDB. A total of 7 B cell epitopes were predicted. Amino acid substitutions L75F, T85A, T94A and T245A were part of these B cell epitopes (Table 2). One experimentally validated T cell epitope (MHC-I) was obtained from IEDB. Filtering of predicted MHC-I epitopes resulted into 20 overlapping peptides that were concatenated to yield 10 unique non-redundant regions. A total of 42 overlapping epitopes were predicted for MHC-II that resulted in 8 discrete regions. Amino acid substitutions T85A, S122P/A, S134R, T245A, L266F, S269P, S270N and I306L were part of predicted T cell epitopes (Table 3).

Table 2

B cell epitopes (the amino acid and the substitution is highlighted in red and the substitution is listed in parentheses)

Description	Start	End
Experimentally validated B cell epitopes from IEDB
QLYKTCKQAGTCPPDIIPKV	17	36
RTGYIPL(F)GTRPPT	69	81
LVEETSFIDAGAP	108	120
STHNYEEIPMDTFIVSTNPNTVTSSTPI	189	216
SGYIPANTTIPF	391	402
YMLRKRRKRLPYFF	453	466
Predicted B cell epitopes
HKRSAKRTKRASATQLYKTC	3	22
GTRPPTATDT(A)LAPVRPPLT(A)VDPV	76	98
TTHNNPTFTDPSVLQPP	156	172
VARLGLYSRTTQQVKVVDPAFVTT(A)P	222	246
LDIVALHRPALTSRRTG	284	300
GAKVHYYYDLS	321	331
FYLHPSYYMLRKRRKR	446	461

Table 3

T cell epitopes (the amino acid and the substitution is highlighted in red and the substitution is listed in parentheses)

Description	Start	End
Experimentally validated T cell epitopes from IEDB
AILDINNTV	144	152
Predicted MHC-I epitope
ILQYGSMGVFF	45	55
TATDT(A)LAPV	81	89
SLVEETSFIDAGAPTS(P/A)V	107	123
AILDINNTV	144	152
AETGGHFTL	175	183
GLYSRTTQQVKVVDPAF	226	242
KLITYDNPAY	248	257
DFLDIVALHR	282	291
RI(L)GNKQTLRTRSGKSIGAKVHYY	305	327
SYYMLRKRRKRLPYFFSDV	451	469
Predicted MHC-II epitope
GSMGVFFGGLGIGTGSGT	49	66
EETSFIDAGAPTS(P/A)VPSIP	110	127
PDVSGFS(R)ITTSTDTPAILDINNTVTTVT	128	156
IPMDTFIVSTNPNTVTSST	196	214
KVVDPAFVTT(A)PTKLITYDNP	236	255
AYEGIDVDNTL(F)YFS(P)S(N)N	256	271
IAPDPDFLDIVALHRPALTSR	277	297
FYLHPSYYMLRKRRKRLPYFFSDVSLAA	446	473

Discussion

HPV-16 is the predominant genotype worldwide, associated with invasive cervical cancers and hence remains as the focus for HPV diagnostic development and vaccine research. L2 protein of HPV is a major component required for virus assembly along with L1 and helps in transport of viral particles to the host cell nucleus. L2 protein also helps in mediating the increased efficiency of formation of VLPs by binding with L1 [31, 32]. The C-terminal amino acid residues 396–439 of HPV11 L2 are shown to form the L1 binding domain [11, 33]. This L2 region of HPV11 corresponds to 412–455 amino acid residues of HPV-16. Seven non-synonymous amino acid substitutions were observed in this region, of which I420T and A443G substitutions were also reported at the corresponding positions in L2 of HPV11 [11]. These substitutions might affect the binding of L2 with respective L1 molecule and thereby impact formation of VLPs. It is also noteworthy that these substitutions (I420T and A443G) were predominantly observed in samples from high-grade cervical disease and may have a role in cervical oncogenesis as well.

The phylogenetic tree constructed from L2 gene sequences showed that majority of the isolates (92.5%) clustered in accordance with the lineages assigned on the basis of L1 as reported previously [12]. However, changes in the sub-lineage clustering pattern were observed, which indicates the differential evolution of L1 and L2.

The L2 residues 1-MRHKRSAKRTKR-12 and 456-RKRRKR-461 constitute the nuclear localization signal (NLS) and the region 296-SRRTGIRYSRIGNKQTLRTRS-316 constitutes the arginine rich nuclear retention signal (NRS). Substitution of arginine residues within NRS lead to reduction in L2 concentration in nucleus even though both the NLS are retained [34]. The only substitution observed in the NRS was I306L that occurs immediately after arginine at 305th position and might impact L2 concentration in nucleus, which needs to be validated experimentally.

Nucleotide variations and the resultant amino acid substitutions can alter the L2 protein properties which can impact the carcinogenic potential. We observed that substitutions T245A, L266F, S378V and S384A were significantly associated with the high‐grade cervical disease. These substitutions need further functional validation to decipher their role as molecular marker(s) of cervical carcinogenesis. Contrary to a previous report we did not find difference in the prevalence of non-synonymous variations between samples from normal and high-grade cervical disease, though the occurrence of multiple non-synonymous variations differed significantly [13].

In the present study, the potential impact of substitutions on antigenicity and immunogenicity was evaluated. The amino acid residues 69–81 and 108–120 of L2 protein are highly conserved and play an important role in inducing neutralizing antibodies [35, 36]. Two monoclonal antibodies (mAbs) KIL2 and MAb6 recognize L2 regions, 73–79 and 65–81 respectively [37, 38]. The L75F substitution is part of these antibody binding sites and may impact the production of neutralizing antibodies. In addition to experimentally validated epitopes reported in IEDB, we predicted, both B and T cell epitopes using computational methods. L75F, the only experimentally validated substitution in B cell epitope was rarely observed in the 148 samples studied. Recently, a number of HPV vaccine constructs based on L2 covering different epitope regions are reported [39-42]. The substitution T85A observed in our study is part of the L2 region of the proposed therapeutic and prophylactic vaccine construct based on L2, E6 and E7 regions [41]. Therefore, T85A substitution might impact vaccine efficacy being common to both B and MHC-I epitopes.

The substitution S122P/A is common to the epitope region recognized by both MHC-I and MHC-II. The substitution T245A is part of both B and T cell epitopes (recognized by MHC-II), while the substitution L266F is part of epitope recognized by MHC-II. These substitutions were observed to be associated with high-grade cervical lesions. Thus, the regions containing substitutions (T245A and L266F) need to be given consideration for disease severity status, immune response based on antigenic diversity, validating existing vaccines and designing novel vaccine candidates.

To conclude, the findings reported in this study would help to understand the impact of L2 protein substitutions on cervical oncogenesis as well as in considerations for design of next-generation L2-based HPV vaccine, subject to experimental validations.