Targeted, Deep Sequencing Reveals Full Methylation Profiles of Multiple HPV Types and Potential Biomarkers for Cervical Cancer Progression

Persistent infection with high-risk human papillomavirus (HPV) types is essential for the development of cervical cancer and its precursor, cervical intraepithelial neoplasia grade 3 (CIN3). However, not all precancerous lesions progress to invasive cancer. In fact, lesions in up to 30% to 40% of patients with CIN2-3 will spontaneously regress (1, 2), whereas progression of CIN3 to invasive cervical cancer (ICC) is estimated to occur in <1.5% of untreated patients within 24 months (2). Accurate diagnosis of CIN3 versus ICC allows tailored treatment to the cervix and preservation of fertility when feasible. Thus, it is surprising that so few studies aim to identify biomarkers that could potentially distinguish precancerous lesions that will progress to invasive cancer from those that will not. Viral-based markers are attractive candidates due to their unique roles in cervical carcinogenesis.

Methylation status of HPV CpG dinucleotide sites (CpG) is one feature that has been examined in preinvasive and invasive cervical lesions. Although viral genomes lack CpG islands similar to those in the human genome, they do have numerous CpGs that can be methylated by host-derived DNA methyltransferases (3), possibly altering expression of viral genes. To date, studies of HPV 16 suggest that cervical disease severity associates with enhanced methylation at CpGs in the L1 gene and upstream regulatory region (URR; refs. 4–6), and limited studies suggest this pattern may also extend to other HPV types (7–11). Furthermore, methylation analysis of the L1 gene in HPV 16 and 18 may be a potential predictor of CIN to ICC progression (8, 12). Although multiple studies support the use of viral methylation status for cervical cancer screening, many CpGs have not been adequately studied, most likely due to limitations of common methylation assays. In fact, few publications evaluate CpG methylation across the entire HPV 16 genome (13). Thus, a need remains for larger scale studies to characterize the complete DNA methylomes of many HPV types across the disease spectrum.

Here, we sought to expand the knowledge base by comprehensive analysis of methylation profiles of many different HPV types in a cohort of largely CIN3 cases and ICC cases. We developed an HPV-capture approach, enriching for up to 63 different HPV types in DNA from precancerous and tumor tissues. Targeted sequencing data for viral methylation in bisulfite (BS)-treated captured libraries were obtained with high-throughput, next-generation sequencing, and CpGs of interest were further validated with pyrosequencing. Our data support previous findings from published studies of HPV DNA methylation, while also providing full methylation profiles of 13 HPV types and revealing new candidate CpGs that may associate with CIN3 progression.

The study protocol was approved by the Institutional Review Board of the Medical College of Wisconsin (Milwaukee, WI). All samples in this case–case comparison study were obtained at the time of diagnosis before any treatment was administered. Subjects in our whole-genome sequencing study had ICC [n = 26 squamous cell carcinoma (SCC); n = 11 non-SCC] or CIN (n = 2 CIN1; n = 2 CIN2; n = 20 CIN3; Table 1). The term non-SCC is defined as all invasive samples that are not of squamous cell histology. Specimens for pyrosequencing validation were made up of 50 ICC (n = 36 SCC; n = 14 non-SCC) and 61 CIN3 (Supplementary Table S1), and overlap between these two methodologies consisted of 21 cases (10 ICC and 10 CIN3). Fresh tissue was biopsied from the primary tumor (ICC) or collected from cervical biopsies (CIN), snap-frozen in liquid nitrogen, and stored at −80°C. Other CIN specimens were obtained cervical swabs. Briefly, samples were collected from the squamocolumnar junction of the cervix using a spatula and cytobrush and transferred to 1.5 mL RNAlater (Ambion) or ThinPrep collection media for storage at −80°C (14). Results were extracted from pathology reports generated at the time of tissue collection. Review of hematoxylin and eosin (H&E) slides was performed by a gynecologic pathologist to confirm diagnosis.

Tumor cell content of sections from frozen cervical tissue samples embedded in optimal cutting temperature medium (OCT) was confirmed to be >70% by H&E staining. Epithelial tissue with CIN histology was excised from stroma in biopsy specimens. Total genomic DNA (gDNA) was extracted from 5 to 10 mg of tissue or cell pellets (0.5 mL cervical swabs) with the Gentra Puregene Tissue Kit (Qiagen). The protocol was modified to include two washes with 1× PBS to dissolve OCT before cell lysis and to exclude RNase A digestion. Formalin-fixed, paraffin-embedded tissue sections (10 μm) were cut and H&E stained for evaluation of tumor content. Any normal tissue present was excised by macrodissection. Sections were deparaffinized with xylene, washed 3 times with xylene, and dehydrated with ethanol. Tissue was then mechanically homogenized using a microfuge tube pestle in the presence of proteinase K (Qiagen) and DNA isolated using the Gentra Puregene System. Quantification of gDNA samples was performed using the Qubit dsDNA BR Assay (Thermo Fisher Scientific) and quality assessed via UV spectroscopy.

For BS-converted libraries, gDNA (3 μg) was sheared to 300 bp fragments using the Covaris S220 System (Covaris; duty factor 10%, 200 cycles/burst, and peak incident power of 140 W for 80 seconds). DNA fragments were end-repaired and adenylated using the SureSelectXT Methyl-Seq Target Enrichment Kit (Agilent Technologies) and purified using the Agencourt AMPure XP Kit (Beckman Coulter Genomics). Library size and concentration were confirmed with an Agilent BioAnalyzer DNA 1000 LabChip (Agilent).

Samples were hybridized with a custom SureSelectXT Library (Agilent) designed to target 63 different HPV types (Supplementary Table S2). Baits were designed against the genomic sequence of each HPV using the eArrary probe design tool (Agilent). For some sequencing runs, xGen Lockdown Probes (Integrated DNA Technologies) targeting exons of 6 cancer-relevant genes were added to the hybridization reaction in an attempt to increase low-sequencing library yields. These probes, however, did not affect subsequent library yield. Probes were 120mers prepared at 4× and 1× tiling coverage for HPV and cancer gene bait, respectively.

Hybrids were captured using Dynabeads MyOne Streptavidin T1 (Thermo Fisher Scientific; Supplementary Fig. S1), and libraries were BS-converted using the EZ-DNA Methylation Kit (Zymo Research) according to the manufacturer's instructions with the following modifications: 10 μg of Yeast tRNA (Life Technologies) was added to each reaction to optimize yield and incubation times modified to 16 cycles of 95°C for 30 seconds and 50°C for 60 minutes. Library concentration was determined via qPCR using the Agilent NGS Library Quantification Kit (Agilent) and were amplified using 15- to 40-cycle PCR with adapter-specific primers. This large range of PCR cycles required to amplify libraries for sequencing may have been due to factors such as DNA input (3 μg), addition of xGEN Lockdown Probes, sample BS conversion, and/or patient variation in HPV viral load. Concentration and size of PCR products was assessed using an Agilent BioAnalyzer DNA High-Sensitivity LabChip (Agilent). Purified PCR products were indexed (6 base pair tag), diluted to 2 nmol/L, and pooled (n = 16) at equimolar ratios for sequencing on the Illumina HiSeq 2000 at the Human and Molecular Genetics Center's Sequencing Facility at the Medical College of Wisconsin. Image analysis and base calling were performed using the standard Illumina pipeline. Zero mismatches were allowed in the index sequence when demultiplexing the data, and all PCR duplicates were removed before analysis.

Adapter sequences were first removed from the output methylation sequence reads using cutadapt (http://code.google.com/p/cutadapt/). Reads with low base quality (<13) were further trimmed and, if less than 25 base pairs, removed by FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Trimmed sequence reads were then analyzed by MethylCoder (15), a software program that generates per-base methylation profiles given a set of reads from BS-treated samples. A genomic short-read nucleotide alignment program (GSNAP; ref. 16) was used for short-read alignment and mapping quality assignment (Phred-scaled probability that the alignment is incorrect). Reference sequences used were human reference (hg19) and sequences of 63 HPV types (Supplementary Table S2; http://pave.niaid.nih.gov; ref. 17). BS conversion rates were estimated from the proportion of unconverted cytosines at non-CpGs. Methylation rate of each CpG was calculated as percent unconverted cytosines in each HPV sample. To test for CpGs containing methylation levels significantly above the BS conversion error rate, we applied a binomial test using a Benjamini–Hochberg corrected P value cutoff of 0.05. Averaging methylation rates of all CpGs in each gene produced HPV methylation levels by gene. CpGs overlapping more than one gene were counted for both genes.

Sequencing reads were aligned to the 63 HPV genomes covered with our custom bait. Read counts aligning to each HPV type (>100) were recorded and used to calculate the proportion of total HPV reads. Only reads uniquely mapped to HPV genomes (mapping quality > 1) were counted. Data from two sequencing runs were subjected to HPV variant analysis to identify potential barcode misallocation leading to cross-contamination among samples, which can occur during sequencing (18). SeqMan Pro (DNASTAR) generated SNP reports for comparison of HPV variants in samples sequenced on the same Illumina lane. If a sample contained identical HPV variants found in another sample ran on the same lane, while also having no unique variants of its own, it was considered contamination and coded as negative for that HPV type. HPV types were assigned to samples if the proportion of mapped reads for that type was ≥15% of the sample's total HPV read count (Supplementary Fig. S2).

Isolated gDNA (5–10 ng) from all samples was BS converted as described above and subjected to PCR using AmpliTaq Gold DNA Polymerase with 0.1 μmol/L each of forward and reverse primer (Supplementary Table S3), with an initial denature step at 95°C for 10 minutes, 50 cycles of 95°C for 30 seconds, 60°C for 30 seconds, 72°C for 1 minute, and final extension step at 72°C for 10 minutes. PCR products were analyzed via gel electrophoresis and pyrosequencing performed using 10 μL of biotinylated PCR product, Streptavidin Sepharose High Performance Beads (GE Healthcare Bio-Sciences), and Pyromark Gold Q96 Reagent Kit (Qiagen). DNA methylation was analyzed on a PyroMark Q96 MD using PyroMark CpG Software (Qiagen).

Global, genic, and CpG-specific methylation differences were examined by two-sample t statistics and Mann–Whitney tests using GraphPad Prism version 6.05 for Windows (GraphPad Software). To examine methylation variation among HPV types, global and genic methylation data were analyzed using one-way ANOVA (F values) and Tukey multiple comparisons post hoc tests. To evaluate the ability of methylation levels to distinguish precancer from cancer, we estimated the sensitivity and specificity using ROC curves. Area under ROC curve was calculated using the caret package in R (http://www.r-project.org/). Two-sided statistical analyses were executed in R (http://www.r-project.org/).

Samples (Table 1) were subjected to HPV-specific enrichment (Supplementary Fig. S1) and BS treatment. On average, our method generated 18.6 million reads/sample, covered 98.5% of viral genomes, and achieved 4,892× depth-of-sequence coverage. The average BS conversion rate was 99.3%, and approximately 69% of BS-converted reads could be uniquely mapped to reference genomes (human and HPV; Supplementary Table S4). Non-HPV sequences were dispersed randomly throughout the human genome, except in sequencing runs that included xGen Lockdown Probes (see Materials and Methods section).

Short sequencing reads uniquely aligned to each HPV type in each sample were used to calculate the proportion of total HPV reads for each HPV type in that sample. As a result, reads from multiple HPV types were identified in each sample. To accurately type samples, we excluded any HPV type with <100 mapped reads and employed a ≥15% cutoff for the percentage of type-specific HPV reads to total HPV reads (Supplementary Fig. S2 and Methods). We identified 13 unique HPV types among the 61 samples (Supplementary Fig. S3) with each sample positive for 1–3 HPV types (>1 types were more common in CIN than ICC). Furthermore, sequencing data showed an 86.9% correspondence with PCR-based genotyping of the same specimens. The lack of complete correspondence was not surprising because some HPV types, especially more uncommon types, are not captured when using degenerate PCR primers or sequence-specific assays. We found 100% correspondence (10/10 samples) in HPV types assigned using our method compared with RNA sequencing data from The Cancer Genome Atlas (19). HPV 16 was the most represented type in our samples, making it impossible for statistically meaningful comparisons in samples harboring other HPV types.

The 63 HPV genomes targeted in this study contained a median of 159 CpGs per type, ranging from 109 to 290. CpGs were numbered according to genomic reference sequences provided by the Papillomavirus Episteme (http://pave.niaid.nih.gov; ref. 17). Methylation rates of individual CpGs in the four most common HPV types from our samples are illustrated in Figs. 1 and 2 (complete methylation data provided in Supplementary Table S5). Averaging methylation rates across all CpGs for all HPV types, we found that global viral methylation in ICC was 5 times higher than in CIN (****P = 2.7 × 10⁻⁸; Fig. 3,A), and this difference was not affected by removal of CIN1-2 and non-SCC samples (****, P = 2.93 × 10⁻⁷; Fig. 3,B). Finally, global HPV methylation of CIN1-2 only samples was significantly lower compared with ICC (P = 0.0042), SCC (P = 0.0032), or non-SCC samples (P = 0.0142), while no significant difference in global methylation was observed comparing SCC with non-SCC (P = 0.7332; data not shown).

Differences in methylation levels between CIN and ICC (Supplementary Fig. S4A) or CIN3 and SCC (Supplementary Fig. S4B) varied substantially among CpGs within genes of the 13 unique viral genomes. Averaging methylation percentages of all CpGs within each gene, we found that the L1, L2, and E5 genes across all HPV types were highly methylated, accounting for >60% of total methylation in ICC or SCC, whereas the URR and E1, E6, and E7 genes were methylated at a considerably lower rate (<20% of total HPV methylation; Supplementary Fig. S4A and S4B). Because HPV 16 was the most common type in our samples, we further examined methylation at single CpGs in all HPV 16–positive samples. Interestingly, we found that methylation levels at individual CpGs in HPV 16 were significantly higher in ICC than CIN (Supplementary Fig. S5A) and in SCC versus CIN3 (Supplementary Fig. S5B).

To assess whether HPV methylation could accurately predict disease status, we used it as a quantitative classifier, constructing ROC curves to determine sensitivity and specificity. The area under an ROC curve (AUC) reflects how well a classifier performs, with a perfect classification yielding an AUC value of up to 1, whereas a random guess would yield an AUC of 0.5. We performed this analysis in two ways: (i) at the gene level (average methylation of all CpGs/gene) and (ii) at the individual CpG level (HPV 16–positive samples only). Analysis of gene methylation in our combined CIN and ICC samples' 13 HPV types showed differing abilities to classify disease status. L1 gene methylation produced the highest AUC value (up to 0.97), suggesting that L1 methylation as a classifier could have 100% sensitivity and more than 95% specificity for detecting ICC (Fig. 4,A; Supplementary Fig. S6A). With removal of CIN1-2 and non-SCC samples, we observed that L2 gene methylation produced a very similar AUC (0.969) to L1 gene methylation (0.965; Fig. 4,B; Supplementary Fig. S6B). Overall, methylation levels of genes producing the top 5 highest AUCs did not change upon removal of these samples (Fig. 4,A and B).

In evaluating the performance of individual CpGs in HPV 16–positive samples, we found a few CpGs that reached AUC values as high as 0.95. Combining results for more than one CpG or gene did not significantly improve the ability to distinguish CIN from ICC. Of the top 5 CpGs ordered by AUC (Fig. 4,C; Supplementary Fig. S6C), 4 were localized within the E1 gene and, to our knowledge, have never before been associated with ICC disease status. When examining the performance of individual CpGs in HPV 16–positive CIN3 and SCC samples, only one of these top 5 CpGs dropped out (3042 replaced by 1558; Fig. 4,D; Supplementary Fig. S6D), whereas the strongest performer (CpG 1870) improved (AUC = 1) with removal of CIN1-2 and non-SCC subjects.

To confirm differential methylation of the 4 CpGs in the E1 gene of HPV 16, we employed pyrosequencing. Primers were designed to amplify these specific CpGs (Supplementary Table S3); however, five distinct primer designs for CpG 1800 failed to generate a PCR product. Thus, we focused on remaining sites (CpG 978, 1870, and 1958) with 1 additional site (CpG 927) also amplified by primers for CpG 978. BS-converted DNA from 61 CIN3 and 50 ICC samples was subjected to pyrosequencing (Supplementary Table S1). Pyrosequencing results confirmed significantly higher methylation all 4 CpGs in ICC versus CIN3 (Fig. 5,A) and SCC versus CIN3 (Fig. 5,B). Importantly, we did not observe any significant difference in the same samples across methodologies (whole-genome vs. pyrosequencing). The average failure rate for generating a PCR product was 10.7% and 23.3% for ICC and CIN3 samples, respectively, and may be explained by lower quality gDNA resulting from BS treatment, lower viral load, or sequence change due to an integration event. Finally, no significant difference in methylation of these CpGs was observed in SCC compared with non-SCC samples (data not shown).

The current study applied a hybridization capture system enriching for up to 63 HPV types in CIN and ICC specimens. Combined with BS treatment, we examined CpG methylation across entire viral genomes. Our data demonstrate that, regardless of HPV type, average global and gene-specific viral methylation is significantly higher in ICC compared with CIN or in SCC compared with CIN3. We also highlight 5 individual CpGs in HPV 16 that can distinguish CIN3 from ICC or CIN3 from SCC. Four of these 5 CpGs are within the E1 gene of HPV 16, have never before been associated with ICC status, and were validated using pyrosequencing. Importantly, the power of the current study lies with distinguishing CIN3 from ICC due to the fact that our “CIN” group consisted of only 4 CIN1-2 cases. Furthermore, our pyrosequencing validation was limited to a cohort of samples comprised of CIN3 cases. To truly make a meaningful comparison of CIN versus ICC, our findings must be validated in a larger cohort of CIN1-2 samples. Nonetheless, the current study provides a strong basis for further exploration of novel CpGs that associate with ICC status and highlight the need for more comprehensive, whole-genome methylation profiling of HPV types across the disease spectrum.

Although some studies have focused on host biomarkers to differentiate CIN from ICC, viral changes may provide disease-specific markers with less variability in interpretation (20, 21). Unfortunately, many studies examining HPV methylation and its association with cervical disease have focused on specific genes or CpGs in the most common HPV types, HPV 16 and 18, leaving rarer HPV types underrepresented in the literature (22). We provide full genome methylation data from 13 different HPV types in CIN and ICC, several not reported before. Our results demonstrate significantly higher global and gene methylation rates in ICC versus CIN or SCC versus CIN3, confirming previous findings of increased methylation in the L1 and URR regions of HPV 16 and 18 in higher grade lesions and cancer (5–9, 13, 23, 24). In accordance with another study (9), we found that, across 13 different HPV types, higher differences in methylation (CIN vs. ICC or CIN3 vs. SCC) occurred in the L1, L2, and E5 genes, while methylation was considerably lower in the URR and E1, E6, and E7 genes. E2-binding sites (E2BS) in HPV 16 and HPV 18 have been shown to be differentially methylated throughout cervical carcinogenesis, leading to downstream effects on E6/E7 oncogene expression (25). Although our HPV 16- and HPV 18–positive samples did not exhibit this same pattern of E2BS methylation, the specifics underlying this discrepancy are unclear and suggest the need for additional work on E2BS methylation and its association with disease status.

Our findings demonstrate that methylation levels at the HPV gene level or at single CpGs can serve as a quantitative classifier for accurately differentiating CIN3 from ICC. Using two methods, we show that individual methylation rates of 5 CpGs within the E1 gene of HPV 16 can distinguish between CIN3 and ICC or SCC. Intriguingly, this is the first example defining disease status using methylation at these CpGs in HPV 16. Examination of CpG methylation in E1 of HPV 16 is exceedingly limited, even in the most comprehensive studies of HPV methylation (4, 7, 9, 13, 22). In fact, to our knowledge, the only evidence of higher methylation at an E1 CpG correlating with cervical disease was restricted to HPV 31 and used only noncancer controls (9). Although we are not sure why E1 CpGs are typically underrepresented in the literature, our results highlight a vital need for more whole-genome methylation analysis. Finally, further studies examining methylation at these specific CpGs in E1 are required to corroborate their potential to assign ICC status and to further define their association with CIN3 cases that are more likely to progress to ICC.

Differential methylation at E1 CpGs as an indicator of cervical disease progression is a novel and provocative concept. The E1 gene is the most conserved of the HPV genes (26, 27) and encodes a protein essential for viral replication. Disruption of E1 function is a critical step in cervical carcinogenesis and can occur via multiple mechanisms. For example, the E1 and E2 genes are most commonly disrupted by integration into the host genome, leading to transcriptional activation of the E6 and E7 oncogenes (28–32) and loss of viral replication control. Also, E1 can be degraded as part of a negative feedback loop with augmented NF-κB signaling (33) and can also induce DNA damage response to facilitate viral genome amplification when its nuclear export is dysfunctional (34). Thus, E1 hypermethylation may represent a redundant mechanism of E1 silencing in cervical disease progression. Alternatively, hypermethylation at E1 body CpGs could lead to alternative splicing of E1 and dysfunctional protein products that aid in driving carcinogenesis. Finally, although there are many functional binding domains throughout E1 (35), the CpGs studied here reside upstream, making it unlikely that their methylation would hinder important protein binding. More work is needed to fully elucidate the functional consequences of increased E1 CpG methylation.

In closing, we would like to acknowledge limitations of the current study. First, sample multiplexing is a cost-effective and routine method for targeting specific genomic regions or when working with small genomes. However, single-indexing protocols may cause barcode misallocation, leading to sample cross-contamination. In our study, 16 samples were pooled and sequenced in a single lane. Although cross-contamination could affect sequence-based HPV genotyping, it is unlikely to have greatly affected our HPV methylation results. The observed methylation differences were typically large with little variance; thus, any contamination would likely result in insignificant changes in our measurements. In addition, pyrosequencing validated specific CpG methylation results from our whole-genome sequencing approach. Second, although 61 samples were analyzed, most were HPV 16 positive, thus limiting our ability to do more comparisons in samples infected with HPV types other than HPV 16. Third, the majority of samples in this study were CIN3 and SCC, limiting the power of other comparisons, although removal of CIN1-2 and non-SCC samples did not significantly affect our main results. Finally, we observed no difference in methylation of off-target human genomic sequences in ICC versus CIN. Yet the current study was designed to enrich for viral DNA; therefore, the human genomic sequences we could analyze may not represent specific methylation differences in the human genome, such as at promoter regions (e.g., tumor suppressor genes), between ICC and CIN.

No potential conflicts of interest were disclosed.

Conception and design: M. Iden, J.S. Rader

Development of methodology: Y.-W. Huang, J.S. Rader

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): S. Fye, Y.-W. Huang, E. Hopp, J.S. Rader

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): P. Liu, M. Iden, S. Fye, Y.-W. Huang, C. Chu, Y. Lu, J.S. Rader

Writing, review, and/or revision of the manuscript: P. Liu, M. Iden, Y.-W. Huang, J.S. Rader

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): M. Iden, J.S. Rader

Study supervision: J.S. Rader

This work was supported in part by start-up funds from Advancing a Healthier Wisconsin Fund (FP00001701 to P. Liu; and FP00001703 to Y. Lu), the Women's Health Research Program at the Medical College of Wisconsin (to J.S. Rader), and a seed grant from the Institute for Translational Medicine of Zhejiang University and National Natural Science Foundation of China (nr. 81372514 and 81572256 to P. Liu; and 81472420 and 31401125 to Y. Lu).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.