Abstract
Steven F. Gameiro1,*, Bart Kolendowski2,*, Ali Zhang1, John W. Barrett3, Anthony C. Nichols3,4,5, Joe Torchia2,4,5 and Joe S. Mymryk1,3,4,5
1Department of Microbiology and Immunology, The University of Western Ontario, London, ON N6A 3K7, Canada
2Department of Biochemistry, The University of Western Ontario, London, ON N6A 3K7, Canada
3Department of Otolaryngology, Head & Neck Surgery, The University of Western Ontario, London, ON N6A 3K7, Canada
4Department of Oncology, The University of Western Ontario, London, ON N6A 3K7, Canada
5London Regional Cancer Program, Lawson Health Research Institute, London, ON N6C 2R5, Canada
*These authors have contributed equally to this work
Correspondence to:
Joe S. Mymryk, email: [email protected]
Keywords: cancer, human papillomavirus, epigenetics, gene expression, methylation
Received: April 13, 2017 Accepted: July 03, 2017 Published: August 03, 2017
ABSTRACT
High-risk human papillomaviruses (HPV) cause cancer at multiple distinct anatomical locations. Regardless of the tissue of origin, most HPV positive (HPV+) cancers show highly upregulated expression of the p16 product of the cyclin-dependent kinase inhibitor 2A (CDKN2A) gene. Paradoxically, HPV+ tumor cells require continuous expression of this tumor suppressor for survival. Thus, restoration of normal p16 regulation has potential therapeutic value against HPV induced cancers. Normally, p16 transcription is tightly controlled at the epigenetic level via polycomb repressive complex-mediated tri-methylation of histone 3 lysine 27 (H3K27me3). Although a mechanism by which HPV induces p16 has been proposed based on tissue culture models, it has not been extensively validated in human tumors. In this study, we used data from over 800 human cervical and head and neck tumors from The Cancer Genome Atlas (TCGA) to test this model. We determined the impact of HPV status on expression from the CDKN2A locus, the adjacent CDKN2B locus, and transcript levels of key epigenetic regulators of these loci. As expected, HPV+ tumors from both anatomical sites exhibited high levels of p16. Furthermore, HPV+ tumors expressed higher levels of KDM6A, which demethylates H3K27me3. CpG methylation of the CDKN2A locus was also consistently altered in HPV+ tumors. This data validates previous tissue culture studies and identifies remarkable similarities between the effects of HPV on gene expression and DNA methylation in both cervical and oral tumors in large human cohorts. Furthermore, these results support a model whereby HPV-mediated dysregulation of CDKN2A transcription requires KDM6A, a potentially druggable target.
INTRODUCTION
Human papillomaviruses (HPV) that infect mucosal tissues are classified as low- or high-risk based on the frequency with which they are associated with cancer. Although both types dysregulate normal cell growth, infection by high-risk HPV is a causative agent for cancers at multiple distinct anatomical subsites [1, 2]. HPV was first recognized as the cause of virtually all cervical carcinomas, and was then subsequently found to contribute to other anogenital cancers as well [2-5]. More recently, HPV infection was recognized as the agent responsible for an epidemic of oropharyngeal cancers [6, 7]. HPV has also been suggested to cause a subset of lung, prostate, bladder, and breast cancers [8-11]. Taken together, HPV infection is thought to be responsible for at least 5% of human cancers worldwide [12].
HPV encodes two main oncogenes, E6 and E7 [13]. Both are constitutively expressed in HPV positive (HPV+) tumors, and suppression of expression of either of these two viral proteins causes HPV dependent tumor cells to senesce and die [1, 14-19]. Both E6 and E7 perform multiple functions in an infected and/or cancerous cell. These two viral oncoproteins function by interacting with multiple key cellular regulatory proteins to dysregulate gene expression and growth [13]. For example, p53 is bound and degraded by E6, while the retinoblastoma protein (Rb) and family members are important targets of E7 [20, 21]. Rb controls exit from the G1 phase of the cell cycle, and the interaction and subsequent degradation of Rb by E7 induces inappropriate cell cycle progression in HPV infected cells [22, 23]. In high-risk HPV infections, viral oncoproteins function to uncouple cell growth and differentiation and contribute to the formation of epithelial dysplasia, which may progress to carcinoma if the infection is not resolved [13].
The HPV E7 oncoproteins alter the regulation of many host cellular genes and much effort has been devoted to understanding the mechanisms by which this occurs [24]. A subset of the effects of E7 on the cellular transcriptome are mediated by epigenetic changes [25]. Specifically, E7 expression has been reported to reduce the global levels of tri-methylated lysine 27 on histone 3 (H3K27me3), a mark of transcriptionally silenced chromatin [26-28]. Paradoxically, several cell culture studies have shown that HPV induces a significant upregulation in the expression level of the enhancer of zeste homolog 2 (EZH2) component of the polycomb repressive complex 2 (PRC2) [28-30], which is the methyltransferase responsible for mono-, di-, and tri-methylation of H3K27 [27, 31]. Both studies also observed an HPV dependent increase in lysine demethylase 6A (KDM6A) or lysine demethylase 6B (KDM6B)—the demethylases that convert H3K27me3 to the di-methylated or mono-methylated forms [28, 30]. Thus, HPV increases expression of the enzymes responsible for creating and removing H3K27me3, complicating our understanding of how the global decrease in H3K27me3 is achieved in an HPV infection.
Mechanistic studies of these HPV dependent changes on epigenetic regulators of gene expression have focused on altered p16 expression, which is present at high levels in nearly all HPV induced tumors [32, 33]. Indeed, p16 was commonly used as a surrogate marker of HPV status during pathological assessment of tumors, but this method to determine HPV status has largely been superseded by the implementation of molecularly based diagnostics [33-35]. p16 and p14 are products of the cyclin-dependent kinase inhibitor 2A (CDKN2A) gene locus and normally function as negative regulators of the cell cycle. p16 induces cell cycle arrest and senescence by inhibiting E2F-mediated transcription. p14 induces cell cycle arrest and apoptosis by facilitating p53 function [36]. Neither of these products inhibit cell cycle progression in HPV induced cancers, as the p53 and Rb pathways have been abrogated by E6 and E7, respectively [20, 21]. However, transcription from the CDKN2A locus is tightly regulated by epigenetic changes, with PRC2 serving a key role in repressing expression by catalyzing tri-methylation of H3K27 across this locus. Subsequently, the polycomb repressive complex 1 (PRC1) is recruited to maintain a transcriptionally repressive chromatin environment that is typically reflected by a local increase in DNA methylation [31]. Thus, this well-characterized locus serves as a good model of epigenetic dysregulation induced by HPV infection [37, 38].
In a detailed series of tissue culture based experiments, evidence was obtained to suggest that the HPV-mediated reduction in H3K27me3 and subsequent transcriptional activation of p16 required E7 to induce expression of the H3K27 specific demethylases KDM6A and KDM6B [28, 39]. Knockdown of either KDM6A or KDM6B in HPV+ cervical carcinoma cell lines not only reduced p16 expression, but also induced cell death [28, 39]. Subsequent experiments showed that E7 expression in cell lines induced an acute dependence on KDM6B for growth. E7 expressing cells were similarly dependent on p16 expression, as treatment with a KDM6 selective small molecule inhibitor was cytotoxic in multiple cervical carcinoma cell lines [28, 39]. Thus, HPV E7 expressing cells are uniquely vulnerable to a targeted agent affecting an epigenetic regulator of cellular gene expression.
While most of this work was done in tissue culture models, similar changes in expression of H3K27me3 regulators were reported in individual cervical carcinomas [28-30]. An analysis of these effects has not been done in any HPV+ oral cavity tumors or on a large scale for any HPV induced tumor site. As such, it remains an open question as to whether these changes in epigenetic regulators are specific for distinct subsites of HPV induced tumors, and how frequently they occur in clinical samples. This information is clearly relevant, considering that cancerous cells expressing E7 appear preferentially sensitive to KDM6 inhibition [28, 39], which could represent a new therapeutic approach for their treatment.
In this study, we used data from over 800 human cervical [40] and head and neck tumors [41] from The Cancer Genome Atlas (TCGA) to determine how high-risk HPV oncogene expression alters expression of H3K27me3 regulators in human tumors. We assessed the impact of HPV status on the transcript levels of key regulators of H3K27 methylation, including EZH2, KDM6A, and KDM6B. We also examined the effect of HPV status on gene expression and DNA methylation status across the CDKN2A genomic region in these two distinct patient cohorts. We found that HPV+ tumors from both anatomical sites exhibited high levels of p16 and p14. Furthermore, HPV+ tumors also expressed higher levels of EZH2 and KDM6A. Analysis of CpG methylation also detected two sites in the CDKN2A locus with consistently altered methylation in HPV+ tumors versus HPV negative (HPV-) samples. Taken together, this data validates conclusions drawn from tissue culture studies and identifies remarkable similarities between the effects of HPV on gene expression and methylation in both cervical and oral tumors in large human cohorts. These results also provide compelling evidence from human tumors that support the current model of HPV dysregulation of H3K27 methylation and cognate gene expression via the upregulation of KDM6A.
RESULTS
Impact of HPV status on CDKN2A and CDKN2B expression in human tumors
The CDKN2A and adjacent cyclin-dependent kinase inhibitor 2B (CDKN2B) loci encode three distinct tumor suppressors: p16 INK4A, p14 ARF, and p15 INK4B. This gene cluster also encodes a long non-coding RNA transcribed in the antisense direction that is referred to as CDKN2B-AS or antisense non-coding RNA in the INK4 locus (ANRIL) (Figure 1). Expression of all three coding genes from this cluster is silenced by polycomb repressive complexes (PRCs) during normal cell growth [31]. While HPV induced cancers commonly express high levels of the p16 product of the CDKN2A locus [32, 33], the expression levels of p14, p15, and the non-coding CDKN2B-AS RNA have not been studied as extensively. We analyzed the TCGA Illumina HiSeq RNA expression data from the head and neck squamous cell carcinoma (HNSC) and cervical carcinoma (CESC) cohorts for expression of all four genes (Figure 2). As expected, HPV+ samples from both the HNSC and CESC cohorts had greatly increased levels of p16 mRNA expression compared to HPV- tumors and normal control tissue. Like p16, a similar upregulation of p14, p15, and CDKN2B-AS was observed in HPV+ cancers with respect to HPV- cancers. These results indicate that PRC-mediated repression is abrogated across the entire CDKN2A and CDKN2B regions of the INK4 locus by HPV infection. The only major difference appears to be the high level of p15 expression present in the normal control samples for the HNSC cohort, which likely reflects some tissue specific, HPV independent regulation.
Figure 1: Organization of the CDKN2A and CDKN2B loci. The cyclin-dependent kinase inhibitor 2A (CDKN2A) locus encodes two functionally unrelated protein products named p16 INK4A and p14 ARF. The adjacent CDKN2B locus encodes p15 INK4B. This gene cluster also encodes a long non-coding RNA transcribed in the antisense direction named CDKN2B-AS or antisense non-coding RNA in the INK4 locus (ANRIL). The position of transcripts (black lines), exons (black vertical bars), coding regions (thick black vertical bars), and orientation (direction of arrowheads) of the major transcripts in this region of chromosome 9 are indicated.
Figure 2: HPV perturbation of expression of the CDKN2A and CDKN2B genes in human head & neck and cervical carcinomas. Normalized RNA-Seq data extracted from the TCGA database for the HNSC (A) and CESC (B) cohorts for HPV+, HPV-, and normal control tissues. Numbers in brackets refer to the number of samples included in each analysis. * p≤0.05, ** p≤0.01, *** p≤0.001, ns- not significant.
Impact of HPV status on expression of PRC2 components and BMI1 in human tumors
The H3K27 methyltransferase activity of PRC2 relies on the catalytic subunit EZH2, as well as the embryonic ectoderm development (EED) and suppressor of zeste 12 (SUZ12) components. Multiple studies have reported EZH2 upregulation in HPV+ cancer cell lines, with similarly higher expression reported in small studies of HPV+ cervical carcinomas [28-30]. Much less is known about the effect of HPV status on the other PRC2 core components. Analysis of TCGA data reveals significantly elevated expression of all three PRC2 core components compared to normal control tissue in both HNSC and CESC (Figure 3). All but SUZ12 in the CESC cohort are similarly elevated as compared to HPV- cancers. Thus, the core components of the PRC2 methyltransferase are all consistently transcribed at high levels in HPV+ human tumors. This is in stark contrast to the reported decreases in global H3K27me3 [28, 30] and the upregulation of the CDKN2A and CDKN2B loci we observe in HPV+ tumors.
Figure 3: HPV perturbation of expression of polycomb components that may regulate CDKN2A transcription in human head & neck and cervical carcinomas. Normalized RNA-Seq data extracted from the TCGA database for the HNSC (A) and CESC (B) cohorts for HPV+, HPV-, and normal control tissues. Numbers in brackets refer to the number of samples included in each analysis. * p≤0.05, ** p≤0.01, *** p≤0.001, ns - not significant, red bracket indicates a comparison that did not achieve significance with a power value <0.8.
Tri-methylation of H3K27 by PRC2 is proposed to recruit PRC1 via its B lymphoma Mo-MLV insertion region 1 (BMI1) component to maintain a transcriptionally repressive state [27]. A previous study reported that BMI1 expression was reduced by transduction of HPV E6 and E7 into human foreskin keratinocytes (HFKs) [30]. A reduction in BMI1 is consistent with a reduction of repression of the CDKN2A locus by PRC1. BMI1 expression in HPV+ CESC tumors was lower than in HPV- or normal control tissues (Figure 3). However, its expression was generally the same or higher in HPV+ HNSC tumors compared to the HPV- or normal control tissues (Figure 3). The lack of a consistent reduction in BMI1 expression in HPV+ tumors suggest that it is not a key factor in HPV activation of CDKN2A and CDKN2B transcription. This is supported by a smaller independent study that found no correlation between p16 and BMI1 expression in head and neck cancers [42].
Impact of HPV status on expression of KDM6A and KDM6B in human tumors
The KDM6A and KDM6B demethylases convert H3K27me3 to the di-methylated or mono-methylated forms [43]. Multiple studies have reported that upregulation of one or both of the KDM6 demethylases occurs in cell lines when HPV E7 is expressed. Furthermore, knockdown of KDM6B or treatment with a small molecule KDM6 inhibitor reduces p16 expression in these models, suggesting a direct causal relationship between KDM6 activity and p16 transcription [28, 39]. Our analysis of the expression of KDM6A reveals that it is significantly upregulated in HPV+ tumors compared to HPV- tumors in both the HNSC and CESC cohorts (Figure 4). In contrast, KDM6B expression is significantly reduced in HPV+ cancers compared to HPV- tumors in both the HNSC and CESC tumor types (Figure 4). Significant changes in KDM6A and KDM6B expression were similarly noted between HPV+ and normal control tissues for the HNSC samples. However, no significant differences were observed for the CESC samples, which could reflect the limited number of normal control samples in this cohort. Overall, the increased expression of KDM6A suggests that it is more likely to play a role in dysregulation of H3K27me3 mediated regulation of transcription in HPV+ cancers.
Figure 4: HPV perturbation of expression of the lysine demethylases KDM6A and KDM6B in human head & neck and cervical carcinomas. Normalized RNA-Seq data extracted from the TCGA database for the HNSC (A) and CESC (B) cohorts for HPV+, HPV-, and normal control tissues. Numbers in brackets refer to the number of samples included in each analysis. * p≤0.05, ** p≤0.01, *** p≤0.001, ns- not significant, red brackets indicate a comparison that did not achieve significance with a power value <0.8.
Assessment of the DNA methylation status of the CDKN2A and CDKN2B loci in human tumors
HPV-mediated changes in epigenetic regulation of transcription of the CDKN2A and CDKN2B gene would be expected to be reflected in methylation status of CpG dinucleotide sequences in or adjacent to these loci. Indeed, a previous report identified altered methylation at 4 CpG loci in this genomic region that correlated with HPV status [30, 38]. We analyzed the Infinium HumanMethylation450 BeadChip array data for the TCGA HNSC and CESC cohorts for consistent alterations in DNA methylation in this region of the genome (Figure 5). A region of DNA upstream of the p14 transcript (probe cg14069088; see Figure 5A for location) was consistently identified to be hypomethylated in HPV+ tumors in both the HNSC and CESC cohorts compared to HPV- tumors or normal control tissue. No differences in methylation were detected in the remainder of the probes across the transcribed regions of these loci in HPV+ versus HPV- or normal control samples. In contrast, a hypermethylated region of DNA was identified within the 3’-end of the CDKN2A locus (probe cg12840719; see Figure 5A for location) in HPV+ tumors regardless of anatomical location. Thus, altered methylation is detected at the CDKN2A locus in HPV+ tumors, which is consistent with an epigenetic mechanism for transcriptional activation of this locus by HPV infection.
Figure 5: HPV perturbation of DNA methylation at the CDKN2A and CDKN2B loci in human head & neck and cervical carcinomas. HPV+ cancers exhibit markedly altered methylation from HPV- tumors or normal control tissues at two regions flanking the CDKN2A locus. (A) Location of the transcripts (black lines), CpG islands and methylation probes used in this study. Normalized methylation data was extracted from the TCGA database for the HNSC (B) and CESC (C) cohorts. The mean methylation levels for HPV+, HPV-, and normal control tissues for all probes located in the vicinity of the CDKN2A and CDKN2B loci are shown. ** p≤ 0.01, *** p≤ 0.001, **** p ≤ 0.0001.
DISCUSSION
HPV infection and subsequent expression of the E6 and E7 oncogenes leads to a large-scale reprogramming of gene expression. This is likely further exacerbated in high-risk HPV+ tumor cells, which typically express higher levels of these powerful viral regulatory proteins compared to levels found in the context of normal infections. Consequently, high-risk HPV transformed cells exhibit many changes in gene expression that contribute to the acquisition and maintenance of cancerous phenotypes unique to viral transformation [13]. In particular, dysregulation of H3K27me3 homeostasis by the E7 proteins from high-risk HPVs contributes to increased expression of loci normally repressed by PRCs, including the CDKN2A locus encoding p16 INK4A [26, 28, 39]. Interestingly, while strong p16 expression is a near universal hallmark of high-risk HPV+ carcinomas of the cervix or oral cavity [32-34], p16 is not expressed in high-risk HPV+ low-grade lesions and variably expressed in higher grade precancerous lesions [44, 45] or lesions caused by low-risk HPV types [32, 46, 47]. Thus, expression of high-risk HPV oncogenes alone is not sufficient to activate high levels of p16 expression, and additional adaptations acquired during tumor development likely contribute to this process.
While cell culture models have allowed detailed analyses of factors involved in virus mediated activation of p16 gene expression from the CDKN2A locus, far less is known about how similar these effects are to actual HPV+ tumors. In this study, our goal was to determine if the changes in H3K27me3 regulators and their target genes observed in cell culture models of HPV transformation are recapitulated in primary HPV+ tumors from different anatomical locations.
Using data from over 800 primary human tumors, we provide evidence that HPV+ tumors from both the head and neck and cervical sites exhibit remarkably similar alterations in expression of the genes encoded by the CDKN2A locus (Figure 2). It is well established that HPV induces expression of p16 INK4A in both oropharyngeal and cervical carcinomas [33, 34]. While overexpression of the p14 ARF gene from this locus has been observed in these cancers [37, 38], this has not been extensively documented. Our observation that p14 expression, like that of p16, is similarly elevated in comparison to HPV- tumors or normal control tissue clearly demonstrates that HPV activates expression of the entire CDKN2A locus. Furthermore, this activation extends into the adjacent CDKN2B locus, which is also regulated by PRCs [31]. Indeed, both the long non-coding RNA, CDKN2B-AS ANRIL, and p15 INK4B are induced in HPV+ cervical carcinomas with respect to HPV- tumors and normal control tissue; similarly, high levels of both CDKN2B encoded genes are expressed in HPV+ head and neck cancers (Figure 2). Existing literature from independent cohorts confirms that p15 is expressed at lower levels in normal cervical tissue compared to cervical carcinomas [48, 49], supporting our conclusions from the TCGA CESC data. In addition, an independent analysis of a subset of the TCGA HNSC cohort used in our analysis observed similarly high, albeit variable, levels of p15 in HPV+ and normal control samples [50]. Taken together, the similarities between the normalized expression of all four gene products from the CDKN2A and CDKN2B gene loci in HPV+ tumors from both the cervix and oropharynx is extraordinary. These results clearly demonstrate the dominant nature of viral reprogramming of gene expression in these cancers with respect to normal control tissue and HPV- cancers.
A simple mechanism to explain the observed increase of CDKN2A and CDKN2B expression in HPV+ cells would be a decreased expression of the proteins necessary for the H3K27 methyltransferase activity of the PRC2. In tissue culture models, expression of the EZH2 methyltransferase is increased by the presence of HPV oncogenes, and no significant change is observed for the EED and SUZ12 core proteins also required for catalytic activity [28, 30]. Similarly, in HPV+ tumors from both anatomical sites, we find that these genes are expressed at levels equivalent or even higher than normal control tissue and HPV- samples (Figure 3). Furthermore, an analysis of ten biopsies of normal cervical tissue, HPV+ high-grade precancerous lesions, and HPV+ cervical carcinomas detected similarly elevated levels of EZH2 protein in the HPV expressing samples [29]. Thus, there is a remarkably close agreement between the finding reported from tissue culture models and what is observed in human tumors, confirming that a reduction in the overall level of PRC2 expression is not likely responsible for CDKN2A and CDKN2B activation.
Silencing of euchromatin by PRC2 is stabilized and maintained by PRC1. Cell culture experiments suggest that expression of BMI1, a core component of PRC1, is reduced in E6/E7 transduced HFKs [30] and BMI1 expression is inversely correlated with HPV titer in cervical carcinomas [51]. Our analysis of primary human cervical carcinomas indicates that BMI1 expression is significantly reduced in HPV+ tumors with respect to HPV- tumors. However, this is not the case in head and neck tumors, where it is modestly upregulated (Figure 3). For this reason, it seems unlikely that the induction of CDKN2A and CDKN2B expression observed in both these cancers is related to a reduction in expression of the BMI1 component of PRC1. This is further supported by a study reporting that BMI1 protein expression was lost in less than half of the 40 p16 positive oropharyngeal carcinoma oral tumors tested [42]. Thus, a loss of BMI1 expression appears unrelated to increased p16 transcription, at least in the context of oropharyngeal tumors.
An alternative method of activating the CDKN2A and CDKN2B loci could be via the removal of the repressive H3K27me3 modification catalyzed by the PRC2 complex. This histone modification is erased by the KDM6A and KDM6B demethylases. Previous work in tissue culture models established that transduction of HPV E7 into HFKs induced KDM6A and/or KDM6B mRNA and protein levels [28, 30]. Mechanistically, induction of KDM6B was independent of the interaction of E7 with Rb. More importantly, siRNA knockdown of KDM6B antagonized p16 induction by E7, suggesting that upregulation of these demethylases is the mechanism by which HPV reduces global H3K27me3 levels and activates p16 expression [28, 39]. Our analysis of primary human tumors indicates that KDM6A expression is significantly increased in HPV+ head and neck and cervical carcinoma cells with respect to HPV- tumors and likely normal control tissue. However, this is not the case for KDM6B, which is significantly reduced (Figure 4). These data fully support the model proposed from tissue culture experiments, in which enhanced expression of KDM6A, but not KDM6B, by HPV is related to the dysregulation of H3K27me3 homeostasis and subsequent activation of CDKN2A and likely CDKN2B transcription in both these cancer types. One small study reported that expression of both KDM6A and KDM6B is substantially reduced in HPV+ mildly dysplastic lesions; however, stepwise increases were observed for both demethylases in precancerous lesions and cervical carcinoma, such that the carcinomas express both demethylases at levels similar to normal control tissues [52]. This pattern of re-expression correlates well with the acquisition of p16 expression observed during cervical carcinoma progression.
Accompanying the global reduction in H3K27me3, HPV+ head and neck cancers and cervical carcinomas exhibit numerous alterations in DNA methylation [53, 54]. Indeed, we recently reported a methylation signature that predicts HPV status in HNSC samples [55]. Our analysis of methylation data for all the probes across the CDKN2A and CDKN2B loci for the TCGA HNSC and CESC cohorts revealed that there were two consistent alterations in DNA methylation relative to methylation levels present in HPV- tumors and normal tissue (Figure 5). A hypomethylated region of DNA was identified on the south shore of the CpG island located at the start of the p14 transcript in HPV+ tumors in both the HNSC and CESC cohorts. Genome-wide studies have clearly shown that methylation in the immediate vicinity of the transcription start site blocks initiation [56]. As such, the greatly reduced methylation in HPV+ samples at this region is fully consistent with enhanced CDKN2A expression. In addition, we also observed a hypermethylated region of DNA on the north shore of the CpG island closest to the third exon of the p14 and p16 transcripts. Hypermethylation at this location has been described in other studies of HPV+ head and neck and cervical carcinomas using alternative probes [38, 57, 58]. These studies have established that p16 and p14 expression in a given tumor sample is strongly correlated with increased methylation at this position in the CDKN2A locus [38, 57]. Hypermethylation of gene-bodies has recently been shown to occur during elongation through a mechanism whereby RNA Polymerase II recruits Su(var)3-9, Enhancer-of-zeste and Trithorax domain containing 2 (SETD2), a histone methyltransferase, which then tri-methylates H3K36. DNA Methyltransferase 3 (DNMT3) recognizes and binds H3K36me3, leading to the methylation of CpGs [59, 60]. Thus, along with promoter hypomethylation, gene-body methylation correlates well with the observed enhanced transcription from the CDKN2A locus, providing further support for an epigenetic mechanism as the basis for HPV induction of expression from this locus.
Taken together, our analysis of over 350 HPV+ head and neck or cervical carcinomas in comparison with over 450 HPV- cancers from their respective anatomical sites demonstrates that HPV consistently activates expression of both the CDKN2A and CDKN2B loci in these cancers. Interestingly, this extends beyond just p16, with most tumors expressing the other gene products from these loci including p14, ANRIL, and p15. In these cancers, HPV oncogene expression similarly dysregulates expression of key regulators of H3K27me3, including members of PRC2 and the demethylases KDM6A and KDM6B. As a likely result of these and other alterations in epigenetic regulation, the levels of methylation in the center of the CDKN2A and CDKN2B loci are significantly decreased in cells expressing high-risk HPV oncogenes.
Importantly, epigenetic reprogramming of gene expression induced by HPV should be reversible, and could represent a unique therapeutic target. Paradoxically, cervical carcinoma cell lines or other tumor cell lines expressing HPV E7 from several high-risk HPVs appear to require continuous p16 expression for survival. Although p16 normally serves as a tumor suppressor, for as yet unknown reasons, it becomes necessary for continued growth in cells expressing HPV E7. As HPV+ cells require KDM6 demethylase activity to maintain p16 expression, they are sensitive to small molecule inhibitors of KDM6 in cell culture models [28, 39]. However, it has not been determined if actual primary HPV+ head and neck or cervical carcinoma cells will respond similarly. Our analysis of the epigenetic regulators that control p16 expression finds many commonalities between the hundreds of tumors making up the TCGA CESC and HNSC cohorts and cell culture models, strongly suggesting that testing of KDM6 inhibitors should be aggressively pursued as a novel therapy of HPV+ cancers.
MATERIALS AND METHODS
RNA expression comparisons and statistical analysis
Level 3 RSEM normalized Illumina HiSeq RNA expression data for the TCGA head and neck cancer (HNSC) and cervical carcinoma (CESC) cohorts was downloaded from the Broad Genome Data Analysis Centers Firehose server (https://gdac.broadinstitute.org/). Although rare HPV+ samples have been identified in the TCGA bladder urothelial carcinoma (BLCA) and colon adenocarcinoma (COAD) cohorts, they were not present in sufficient numbers for any useful comparisons. For all genes except p16 and p14, the gene level Firehose dataset was used. For p16 and p14, the gene isoform level Firehose dataset was used (uc003zpk.2 and uc003zpl.2 respectively) to adequately discriminate between these two different products of the CDKN2A gene. Normalized expression data was extracted into Microsoft Excel and the HPV status was manually curated based on published datasets [40, 41, 61, 62]. For each gene analyzed, primary patient samples with known HPV status were grouped as HPV+, HPV-, or normal control tissue. Patient samples with unknown HPV status were omitted from our calculations, as were samples obtained from secondary metastatic lesions. This resulted in 73 HPV+, 442 HPV-, and 43 normal control samples with data available for the HNSC gene expression analysis and 278 HPV+, 19 HPV-, and 3 normal control samples available for the cervical carcinoma gene expression analysis. Boxplot comparison of gene expression was performed using BoxPlotR (http://shiny.chemgrid.org/boxplotr/) and assembled into final form using CorelDRAW. For the box plots, center lines show the medians, box limits indicate the 25th and 75th percentiles as determined by R software and whiskers extend 1.5 times the interquartile range from the 25th and 75th percentiles. Statistical significance was calculated using Graphpad Prism v6.01. p-values were assigned using a one-tailed non-parametric Mann-Whitney U test. Post-hoc power calculations were performed with G*Power software version 3.1.9.2 [63], using post-hoc t-test family calculations, with effect size selected as 0.8 and α = 0.05. All comparisons achieved a power value >0.8, or demonstrated significant differences, unless otherwise noted in the text.
DNA methylation comparisons and statistical analysis
Level 3 Infinium HumanMethylation450 BeadChip array data for the TCGA HNSC and CESC cohorts was downloaded from the Broad Genome Data Analysis Centers Firehose server (https://gdac.broadinstitute.org/). 14 methylation probes were identified to be present in the CDKN2A and adjacent CDKN2B region: cg00718440, cg02400248, cg03079681, cg04026675, cg07562918, cg08390209, cg10848754, cg12840719, cg13601799, cg13926295, cg14069088, cg14430974, cg19133618, cg19481686. Methylation data for all patient samples for these probes were extracted from the Firehose files and assembled into a table using Microsoft Excel. HPV status was manually curated based on published datasets [40, 41, 61, 62]. For each gene analyzed, primary patient samples with known HPV status were grouped as HPV+, HPV-, or normal control tissue. Patient samples with unknown HPV status were omitted from our calculations, as were samples obtained from secondary metastatic lesions. This resulted in 73 HPV+ tumors, 442 HPV- tumors, and 43 normal control samples with data available for the HNSC methylation analysis and 278, 19, and 3 samples respectively available for the CESC methylation analysis. The average methylation Beta-value and standard deviation was calculated for each probe for each sample type in each cohort. Data was plotted using Excel and assembled into final form using CorelDRAW. For statistical analysis, Beta-values were converted to M-values using the following equation to improve homoscedasticity of the data:
Author contributions
SFG participated in study design and analyzed the data. BK participated in study design and analyzed the data. AZ analyzed the data and helped prepare figures. JWB, ACN, and JT participated in study design. JSM conceived the study and was responsible for the final interpretation and drafting of the manuscript. All authors contributed to the drafting of the manuscript.
ACKNOWLEDGMENTS
We thank Drs. Scott Bratman, James Pipas, and Paul Cantalupo for providing information regarding HPV status of various TCGA samples and Charles Ishak for helpful discussion.
CONFLICTS OF INTEREST
Authors declare no conflicts of interest.
FUNDING
This work was supported in part from grants provided by the Canadian Institutes of Health Research (IG1-1345172 and MOP#142491) to JSM and ACN. BK is supported by Studentships from the Translational Breast Cancer Research Unit of the London Regional Cancer Program, funded by the Breast Cancer Society of Canada.
REFERENCES
1. zur Hausen H. Papillomavirus infections--a major cause of human cancers. Biochim Biophys Acta. 1996; 1288: F55–78.
2. zur Hausen H. Papillomaviruses and cancer: from basic studies to clinical application. Nat Rev Cancer. 2002; 2: 342–50. https://doi.org/10.1038/nrc798.
3. Crosbie EJ, Einstein MH, Franceschi S, Kitchener HC. Human papillomavirus and cervical cancer. Lancet. 2013; 382: 889–99. https://doi.org/10.1016/S0140-6736(13)60022-7.
4. Guimerà N, Alemany L, Halec G, Pawlita M, Wain GV, Vailén JS, Azike JE, Jenkins D, de Sanjosé S, Quint W, Bosch FX. Human papillomavirus 16 is an aetiological factor of scrotal cancer. Br J Cancer. 2017; 146: 1299. https://doi.org/10.1038/bjc.2017.74.
5. Alemany L, Cubilla A, Halec G, Kasamatsu E, Quirós B, Masferrer E, Tous S, Lloveras B, Hernández-Suarez G, Lonsdale R, Tinoco L, Alejo M, Alvarado-Cabrero I, et al. Role of human papillomavirus in penile carcinomas worldwide. Eur Urol. 2016; 69: 953–61. https://doi.org/10.1016/j.eururo.2015.12.007.
6. Gillison ML, Koch WM, Capone RB, Spafford M, Westra WH, Wu L, Zahurak ML, Daniel RW, Viglione M, Symer DE, Shah KV, Sidransky D. Evidence for a causal association between human papillomavirus and a subset of head and neck cancers. J Natl Cancer Inst. 2000; 92: 709–20. https://doi.org/10.1093/jnci/92.9.709.
7. Syrjanen S. Human papillomavirus (HPV) in head and neck cancer. J Clin Virol. 2005; 32: S59–66. https://doi.org/10.1016/j.jcv.2004.11.017.
8. Singh N, Hussain S, Kakkar N, Singh SK, Sobti RC, Bharadwaj M. Implication of high risk human papillomavirus HR-HPV infection in prostate cancer in Indian population--a pioneering case-control analysis. Sci Rep. 2015; 5: 7822. https://doi.org/10.1038/srep07822.
9. The Cancer Genome Atlas Research Network. Comprehensive molecular characterization of urothelial bladder carcinoma. Nature. 2014; 507: 315–22. https://doi.org/10.1038/nature12965.
10. Lawson JS, Glenn WK, Salyakina D, Delprado W, Clay R, Antonsson A, Heng B, Miyauchi S, Tran DD, Ngan CC, Lutze-Mann L, Whitaker NJ. Human papilloma viruses and breast cancer. Front Oncol. 2015; 5: 515. https://doi.org/10.3389/fonc.2015.00277.
11. Zhai K, Ding J, Shi HZ. HPV and lung cancer risk: a meta-analysis. J Clin Virol. 2015; 63: 84–90. https://doi.org/10.1016/j.jcv.2014.09.014.
12. Forman D, de Martel C, Lacey CJ, Soerjomataram I, Lortet-Tieulent J, Bruni L, Vignat J, Ferlay J, Bray F, Plummer M, Franceschi S. Global burden of human papillomavirus and related diseases. Vaccine. 2012; 30: F12–23. https://doi.org/10.1016/j.vaccine.2012.07.055.
13. Doorbar J, Quint W, Banks L, Bravo IG, Stoler M, Broker TR, Stanley MA. The biology and life-cycle of human papillomaviruses. Vaccine. 2012; 30: F55–70. https://doi.org/10.1016/j.vaccine.2012.06.083.
14. Deng Z, Hasegawa M, Kiyuna A, Matayoshi S, Uehara T, Agena S, Yamashita Y, Ogawa K, Maeda H, Suzuki M. Viral load, physical status, and E6/E7 mRNA expression of human papillomavirus in head and neck squamous cell carcinoma. Head Neck. 2013; 35: 800–8. https://doi.org/10.1002/hed.23034.
15. Goodwin EC, Yang E, Lee CJ, Lee HW, DiMaio D, Hwang ES. Rapid induction of senescence in human cervical carcinoma cells. Proc Natl Acad Sci U S A. 2000; 97: 10978–83. https://doi.org/10.1073/pnas.97.20.10978.
16. Goodwin EC, DiMaio D. Repression of human papillomavirus oncogenes in HeLa cervical carcinoma cells causes the orderly reactivation of dormant tumor suppressor pathways. Proc Natl Acad Sci U S A. 2000; 97: 12513–8. https://doi.org/10.1073/pnas.97.23.12513.
17. Rampias T, Sasaki C, Weinberger P, Psyrri A. E6 and E7 gene silencing and transformed phenotype of human papillomavirus 16-positive oropharyngeal cancer cells. J Natl Cancer Inst. 2009; 101: 412–23. https://doi.org/10.1093/jnci/djp017.
18. Sima N, Wang W, Kong D, Deng D, Xu Q, Zhou J, Xu G, Meng L, Lu Y, Wang S, Ma D. RNA interference against HPV16 E7 oncogene leads to viral E6 and E7 suppression in cervical cancer cells and apoptosis via upregulation of Rb and p53. Apoptosis. 2008; 13: 273–81. https://doi.org/10.1007/s10495-007-0163-8.
19. Jiang M, Milner J. Selective silencing of viral gene expression in HPV-positive human cervical carcinoma cells treated with siRNA, a primer of RNA interference. Oncogene. 2002; 21: 6041–8. https://doi.org/10.1038/sj.onc.1205878.
20. Dyson N, Guida P, Munger K, Harlow E. Homologous sequences in adenovirus E1A and human papillomavirus E7 proteins mediate interaction with the same set of cellular proteins. J Virol. 1992; 66: 6893–902.
21. Kessis TD, Slebos RJ, Nelson WG, Kastan MB, Plunkett BS, Han SM, Lorincz AT, Hedrick L, Cho KR. Human papillomavirus 16 E6 expression disrupts the p53-mediated cellular response to DNA damage. Proc Natl Acad Sci U S A. 1993; 90: 3988–92. https://doi.org/10.1073/pnas.90.9.3988.
22. Brehm A, Miska EA, McCance DJ, Reid JL, Bannister AJ, Kouzarides T. Retinoblastoma protein recruits histone deacetylase to repress transcription. Nature. 1998; 391: 597–601. https://doi.org/10.1038/35404.
23. Boyer SN, Wazer DE, Band V. E7 protein of human papilloma virus-16 induces degradation of retinoblastoma protein through the ubiquitin-proteasome pathway. Cancer Res. 1996; 56: 4620–4.
24. Munger K, Basile JR, Duensing S, Eichten A, Gonzalez SL, Grace M, Zacny VL. Biological activities and molecular targets of the human papillomavirus E7 oncoprotein. Oncogene. 2001; 20: 7888–98. https://doi.org/10.1038/sj.onc.1204860.
25. Durzynska J, Lesniewicz K, Poreba E. Human papillomaviruses in epigenetic regulations. Mutat Res Rev Mutat. 2017; 772: 36–50. https://doi.org/10.1016/j.mrrev.2016.09.006.
26. McLaughlin-Drubin ME, Huh KW, Münger K. Human papillomavirus type 16 E7 oncoprotein associates with E2F6. J Virol. 2008; 82: 8695–705. https://doi.org/10.1128/JVI.00579-08.
27. Schwartz YB, Pirrotta V. Polycomb silencing mechanisms and the management of genomic programmes. Nat Rev Genet. 2007; 8: 9–22. https://doi.org/10.1038/nrg1981.
28. McLaughlin-Drubin ME, Crum CP, Münger K. Human papillomavirus E7 oncoprotein induces KDM6A and KDM6B histone demethylase expression and causes epigenetic reprogramming. Proc Natl Acad Sci U S A. 2011; 108: 2130–5. https://doi.org/10.1073/pnas.1009933108.
29. Holland D, Hoppe-Seyler K, Schuller B, Lohrey C, Maroldt J, Dürst M, Hoppe-Seyler F. Activation of the enhancer of zeste homologue 2 gene by the human papillomavirus E7 oncoprotein. Cancer Res. 2008; 68: 9964–72. https://doi.org/10.1158/0008-5472.CAN-08-1134.
30. Hyland PL, McDade SS, McCloskey R, Dickson GJ, Arthur K, McCance DJ, Patel D. Evidence for alteration of EZH2, BMI1, and KDM6A and epigenetic reprogramming in human papillomavirus type 16 E6/E7-expressing keratinocytes. J Virol. 2011; 85: 10999–1006. https://doi.org/10.1128/JVI.00160-11.
31. Popov N, Gil J. Epigenetic regulation of the INK4b-ARF-INK4a locus: in sickness and in health. Epigenetics. 2010; 5: 685–90. https://doi.org/10.4161/epi.5.8.12996.
32. Sano T, Oyama T, Kashiwabara K, Fukuda T, Nakajima T. Expression status of p16 protein is associated with human papillomavirus oncogenic potential in cervical and genital lesions. Am J Pathol. 1998; 153: 1741–8. https://doi.org/10.1016/S0002-9440(10)65689-1.
33. Klaes R, Friedrich T, Spitkovsky D, Ridder R, Rudy W, Petry U, Dallenbach-Hellweg G, Schmidt D, Knebel Doeberitz von M. Overexpression of p16(INK4A) as a specific marker for dysplastic and neoplastic epithelial cells of the cervix uteri. Int J Cancer. 2001; 92: 276–84. https://doi.org/10.1002/ijc.1174.
34. Hoffmann M, Ihloff AS, Gorogh T, Weise JB, Fazel A, Krams M, Rittgen W, Schwarz E, Kahn T. p16(INK4a) overexpression predicts translational active human papillomavirus infection in tonsillar cancer. Int J Cancer. 2010; 127: 1595–602. https://doi.org/10.1002/ijc.25174.
35. Burd EM. Human papillomavirus laboratory testing: the changing paradigm. Clin Microbiol Rev. 2016; 29: 291–319. https://doi.org/10.1128/CMR.00013-15.
36. Zhang Y, Xiong Y, Yarbrough WG. ARF promotes MDM2 degradation and stabilizes p53: ARF-INK4a locus deletion impairs both the Rb and p53 tumor suppression pathways. Cell. 1998; 92: 725–34. https://doi.org/10.1016/S0092-8674(00)81401-4.
37. Kanao H, Enomoto T, Ueda Y, Fujita M, Nakashima R, Ueno Y, Miyatake T, Yoshizaki T, Buzard GS, Kimura T, Yoshino K, Murata Y. Correlation between p14ARF/p16INK4A expression and HPV infection in uterine cervical cancer. Cancer Lett. 2004; 213: 31–7. https://doi.org/10.1016/j.canlet.2004.03.030.
38. Schlecht NF, Ben-Dayan M, Anayannis N, Lleras RA, Thomas C, Wang Y, Smith RV, Burk RD, Harris TM, Childs G, Ow TJ, Prystowsky MB, Belbin TJ. Epigenetic changes in the CDKN2A locus are associated with differential expression of P16INK4A and P14ARF in HPV-positive oropharyngeal squamous cell carcinoma. Cancer Med. 2015; 4: 342–53. https://doi.org/10.1002/cam4.374.
39. McLaughlin-Drubin ME, Park D, Münger K. Tumor suppressor p16INK4A is necessary for survival of cervical carcinoma cell lines. Proc Natl Acad Sci U S A. 2013; 110: 16175–80. https://doi.org/10.1073/pnas.1310432110.
40. The Cancer Genome Atlas Research Network. Integrated genomic and molecular characterization of cervical cancer. Nature. 2017; 543: 378–84. https://doi.org/10.1038/nature21386.
41. The Cancer Genome Atlas Research Network. Comprehensive genomic characterization of head and neck squamous cell carcinomas. Nature. 2015; 517: 576–82. https://doi.org/10.1038/nature14129.
42. Lundberg M, Renkonen S, Haglund C, Mattila PS, Leivo I, Hagström J, Mäkitie AA. Association of BMI-1 and p16 as prognostic factors for head and neck carcinomas. Acta Otolaryngol. 2016; 136: 501–5. https://doi.org/10.3109/00016489.2015.1122227.
43. Swigut T, Wysocka J. H3K27 demethylases, at long last. Cell. 2007; 131: 29–32. https://doi.org/10.1016/j.cell.2007.09.026.
44. Mills AM, Paquette C, Castle PE, Stoler MH. Risk stratification by p16 immunostaining of CIN1 biopsies: a retrospective study of patients from the quadrivalent HPV vaccine trials. Am J Surg Pathol. 2015; 39: 611–7. https://doi.org/10.1097/PAS.0000000000000374.
45. Woo SB, Cashman EC, Lerman MA. Human papillomavirus-associated oral intraepithelial neoplasia. Mod Pathol. 2013; 26: 1288–97. https://doi.org/10.1038/modpathol.2013.70.
46. Geißler C, Tahtali A, Diensthuber M, Gassner D, Stöver T, Wagenblast J. The role of p16 expression as a predictive marker in HPV-positive oral SCCHN--a retrospective single-center study. Anticancer Res. 2013; 33: 913–6.
47. Mooren JJ, Gültekin SE, Straetmans JM, Haesevoets A, Peutz-Kootstra CJ, Huebbers CU, Dienes HP, Wieland U, Ramaekers FC, Kremer B, Speel EJ, Klussmann JP. P16(INK4A) immunostaining is a strong indicator for high-risk-HPV-associated oropharyngeal carcinomas and dysplasias, but is unreliable to predict low-risk-HPV-infection in head and neck papillomas and laryngeal dysplasias. Int J Cancer. 2014; 134: 2108–17. https://doi.org/10.1002/ijc.28534.
48. Feng W, Xiao J, Zhang Z, Rosen DG, Brown RE, Liu J, Duan X. Senescence and apoptosis in carcinogenesis of cervical squamous carcinoma. Mod Pathol. 2007; 20: 961–6. https://doi.org/10.1038/modpathol.3800927.
49. Zhang Y, Guo L, Xing P, Chen Y, Li F, Zhu W, Lu X. Increased expression of oncogene-induced senescence markers during cervical squamous cell cancer development. Int J Clin Exp Pathol. 2014; 7: 8911–6.
50. Salyakina D, Tsinoremas NF. Non-coding RNAs profiling in head and neck cancers. NPJ Genom Med. 2016; 1: 38. https://doi.org/10.1038/npjgenmed.2015.4.
51. Kim BW, Cho H, Ylaya K, Chung JY, Hewitt SM. Over-expression of BMI1 is associated with favorable prognosis in cervical cancer. Int J Clin Exp Pathol. 2016; 9: 1849–57.
52. Iancu IV, Botezatu A, Plesa A, Huica I. Histone lysine demethylases as epigenetic modifiers in HPV-induced cervical neoplasia. Rom Biotech Lett. 2015; 20: 10236–44.
53. van Kempen PM, Noorlag R, Braunius WW, Stegeman I, Willems SM, Grolman W. Differences in methylation profiles between HPV-positive and HPV-negative oropharynx squamous cell carcinoma: a systematic review. Epigenetics. 2014; 9: 194–203. https://doi.org/10.4161/epi.26881.
54. Lorincz AT. Virtues and weaknesses of DNA methylation as a test for cervical cancer prevention. Acta Cytol. 2016; 60: 501–12. https://doi.org/10.1159/000450595.
55. Papillon-Cavanagh S, Lu C, Gayden T, Mikael LG, Bechet D, Karamboulas C, Ailles L, Karamchandani J, Marchione DM, Garcia BA, Weinreb I, Goldstein D, Lewis PW, et al. Impaired H3K36 methylation defines a subset of head and neck squamous cell carcinomas. Nat Genet. 2017; 49: 180–5. https://doi.org/10.1038/ng.3757.
56. Jones PA. Functions of DNA methylation: islands, start sites, gene bodies and beyond. Nat Rev Genet. 2012; 13: 484–92. https://doi.org/10.1038/nrg3230.
57. Ben-Dayan MM, Ow TJ, Belbin TJ, Wetzler J, Smith RV, Childs G, Diergaarde B, Hayes DN, Grandis JR, Prystowsky MB, Schlecht NF. Nonpromoter methylation of the CDKN2A gene with active transcription is associated with improved locoregional control in laryngeal squamous cell carcinoma. Cancer Med. 2017; 6: 397–407. https://doi.org/10.1002/cam4.961.
58. Wijetunga NA, Belbin TJ, Burk RD, Whitney K, Abadi M, Greally JM, Einstein MH, Schlecht NF. Novel epigenetic changes in CDKN2A are associated with progression of cervical intraepithelial neoplasia. Gynecol Oncol. 2016; 142: 566–73. https://doi.org/10.1016/j.ygyno.2016.07.006.
59. Neri F, Rapelli S, Krepelova A, Incarnato D, Parlato C, Basile G, Maldotti M, Anselmi F, Oliviero S. Intragenic DNA methylation prevents spurious transcription initiation. Nature. 2017; 543: 72–7. https://doi.org/10.1038/nature21373.
60. Baubec T, Colombo DF, Wirbelauer C, Schmidt J, Burger L, Krebs AR, Akalin A, Schübeler D. Genomic profiling of DNA methyltransferases reveals a role for DNMT3B in genic methylation. Nature. 2015; 520: 243–7. https://doi.org/10.1038/nature14176.
61. Bratman SV, Bruce JP, O’Sullivan B, Pugh TJ, Xu W, Yip KW, Liu FF. Human papillomavirus genotype association with survival in head and neck squamous cell carcinoma. JAMA Oncol. 2016; 2: 823–6. https://doi.org/10.1001/jamaoncol.2015.6587.
62. Banister CE, Liu C, Pirisi L, Creek KE, Buckhaults PJ. Identification and characterization of HPV-independent cervical cancers. Oncotarget. 2017; 8: 13375–86. https://doi.org/10.18632/oncotarget.14533.
63. Faul F, Erdfelder E, Lang AG, Buchner A. G*Power 3: a flexible statistical power analysis program for the social, behavioral, and biomedical sciences. Behav Res Methods. 2007; 39: 175–91.
64. Du P, Zhang X, Huang CC, Jafari N, Kibbe WA, Hou L, Lin SM. Comparison of Beta-value and M-value methods for quantifying methylation levels by microarray analysis. BMC Bioinformatics. 2010; 11: 587. https://doi.org/10.1186/1471-2105-11-587.
