Abstract
Genetic abnormalities of cholangiocarcinoma have been widely studied; however, epigenomic changes related to cholangiocarcinogenesis have been less well characterized. We have profiled the DNA methylomes of 28 primary cholangiocarcinoma and six matched adjacent normal tissues using Infinium's HumanMethylation27 BeadChips with the aim of identifying gene sets aberrantly and epigenetically regulated in this tumor type. Using a linear model for microarray data, we identified 1610 differentially methylated autosomal CpG sites, with 809 hypermethylated (representing 603 genes) and 801 hypomethylated (representing 712 genes) in cholangiocarcinoma versus adjacent normal tissues (false-discovery rate ≤ 0.05). Gene ontology and gene set enrichment analyses identified gene sets significantly associated with hypermethylation at linked CpG sites in cholangiocarcinoma including homeobox genes and target genes of PRC2, EED, SUZ12, and histone H3 trimethylation at lysine 27. We confirmed frequent hypermethylation at the homeobox genes HOXA9 and HOXD9 by bisulfite pyrosequencing in a larger cohort of cholangiocarcinoma (n = 102). Our findings indicate a key role for hypermethylation of multiple CpG sites at genes associated with a stem cell-like phenotype as a common molecular aberration in cholangiocarcinoma. These data have implications for cholangiocarcinogenesis, as well as possible novel treatment options using histone methyltransferase inhibitors. Cancer Prev Res; 6(12); 1348–55. ©2013 AACR.
Introduction
Cholangiocarcinoma is a hepatic malignancy originating from the biliary tract, which accounts for 10% to 20% of primary liver cancers (1, 2). Depending on the anatomical site of the tumor, cholangiocarcinomas are classified into being either of extrahepatic or intrahepatic origin (3). The global incidence rates of intrahepatic cholangiocarcinoma, in particular, have markedly increased worldwide (3, 4). The geographic distribution of risk factors account for geographic differences in incidence rates, with the highest incidence reported in the endemic area of liver-fluke infection in northeastern Thailand. In this region, the age-standardized incidence rate (ASIR) is up to 317.6 per 100,000 people, which exceeds the ASIR of less than 2 per 100,000 people in Western countries by more than a hundred-fold (5–9). However, incidences of cholangiocarcinoma have also been steadily increasing in Western countries such as the United Kingdom or United States over the last 40 years (3). Cholangiocarcinomas are particularly difficult to diagnose, and often present at a late clinical stage, which drastically reduces the chances for successful treatment. Surgical resection has been reported to be curative; however, only very few patients are eligible for this treatment option (10). The vast majority of patients will receive nonsurgical therapeutic regimens such as chemotherapy, with low response rates observed overall.
Epigenetic changes in tumors are associated with cancer development and progression (11). Aberrant DNA methylation occurs very early during carcinogenesis before the onset of malignancy, and is thought to interact with genetic aberrations in driving the tumorigenic phenotype of cancer cells (12). Epigenetic silencing at CpG islands, in particular, is observed in many solid tumor types. Aberrant methylation of specific candidate genes in cholangiocarcinoma has been described by several studies (13–16); however, there have been no studies examining methylation profiles throughout the genome. We have quantified DNA methylation levels at a genome-wide scale using Illumina's Infinium Human Methylation27 BeadChip. The Methylation BeadChip covers 27,578 CpG sites across 14,485 consensus coding sequences (CCDS) and well-known cancer genes in the human genome (17).
It has been proposed that there are subpopulations of cells present within tumors that are responsible for sustaining their growth and are phenotypically distinct from the majority of tumor cells (18). These populations of tumor-sustaining cells have characteristics of stem cells, which include prolonged proliferative capacity in vitro and in vivo, expression of specific subsets of genes, and resistance to DNA damage (19, 20). Tumors often lack the terminal differentiation traits possessed by their normal cell states. These parallels have supported the hypothesis that cancer cells arise from undifferentiated stem or progenitor cells or, alternatively, that cancer cells can undergo progressive dedifferentiation during their development (21–23). Furthermore, it has been suggested that genes involved in regulation of a stem cell state may be more vulnerable to aberrant DNA methylation (24). We have, therefore, addressed whether there is enrichment for aberrant DNA methylation of genes associated with a stem cell-like state in cholangiocarcinoma.
Materials and Methods
Tumor samples and DNA extraction
Twenty-eight primary intrahepatic cholangiocarcinoma and six matched adjacent normal samples were supplied by the Liver Fluke and Cholangiocarcinoma Research Center, Khon Kaen University (Khon Kaen, Thailand). Ethical approval was obtained for use of all samples (Ethics Committee of Khon Kaen University - HE510651). The histopathologic features of the tissues are summarized in Supplementary Table S1. Genomic DNA was extracted using the DNeasy Tissue Kit (Qiagen) according to the manufacturer's recommendations.
DNA methylation profiling
Array-based DNA methylation profiling was performed using Infinium Human Methylation27 BeadChips (Illumina), which simultaneously interrogate 27,578 CpG sites within the promoter regions of more than 14,000 genes (15,369 CpG sites located within CpG islands, 12,209 CpG sites located outside of CpG islands, as defined by Gardiner-Garden and Frommer; ref. 25). One microgram genomic DNA was bisulfite modified using the EZ DNA Methylation Kit (Zymo Research, Cambridge, UK). Then, 200 ng of converted DNA were further processed to run BeadArrays according to the manufacturer's instructions. The arrays were imaged using a BeadArray Reader. Each locus is represented by fluorescent signals from two bead types corresponding to the methylated (M) and unmethylated (U) alleles, respectively. The raw signals of unmethylated and methylated bead types (average, 15–30 beads of each type per locus) were background corrected and computed into a β value using the BeadStudio software Methylation Module version 1.0.5 (Illumina). The β value represents the ratio of the intensity of the methylated bead type to the combined locus intensity and reflects the methylation status of a specific CpG site. Beta values are continuous variables ranging from 0 (unmethylated) to 1 (fully methylated), and these were quartile normalized and a subsequent analysis was carried out in R version 2.10.1.
The robustness and reproducibility of the Infinium HumanMethylation27 BeadChips were verified by using five technical replicates of the ovarian cancer cell line PEO1 for correlation experiments. Beta values of these replicates were highly correlated with an average R2 of 0.989 (range, 0.988–0.990; P < 0.00001).
Bisulfite modification and pyrosequencing analysis
Validation of candidate loci at the two HOX genes, HOXA9 and HOXD9, was performed in a bigger panel of Thai cholangiocarcinoma (n = 102) and adjacent normal (n = 24) samples using bisulfite pyrosequencing. One microgram of genomic DNA was bisulfite modified using the EpiTect Bisulfite Kit (Qiagen) according to the manufacturer's instructions. Up to 100 ng of modified DNA were used as a template in subsequent pyrosequencing reactions. Pyrosequencing primer sets covering differentially methylated CpG sites in the HOXA9 and HOXD9 gene promoters were HOXA9_PYRO_F: 5′-Biotin-GGTGATGGTTATTATTGGGGTTT-3′, HOXA9_PYRO_R: 5′- AAAAACTAACCCAAAATCCC-3′, HOXA9_PYRO_S: 5′-CTAACCCAAAATCCCC-3′, HOXD9_PYRO_F: 5′-GGGTATATGGATTTGGGTTTG-3′, HOXD9_PYRO_R: 5′-Biotin-TCCCACCCAACATTACACTATC-3′, HOXD9_PYRO_S: 5′-TAATTATAGTAGTTTTTAGTTTA-3′. Pyrosequencing PCRs were performed in duplicate for each sample in a 25-μL volume containing 0.2 μL FastStart Taq DNA polymerase (Roche Applied Science), 2.5 μL FastStart Buffer (Roche Applied Science), 0.2 mmol/L dNTPs (Applied Biosystems), 0.2 μmol/L primers (each), and 2 μL of modified DNA template using the following conditions: 95°C for 6 minutes, 35 cycles at 95°C for 30 seconds, 63°C (HOXA9) or 60°C (HOXD9) for 30 seconds and 72°C for 30 seconds followed by final extension at 72°C for 5 minutes. Methylation frequencies were determined by pyrosequencing of PCR products using the PSQ HS96 Pyrosequencing System (Biotage) and converting the resulting pyrograms into numerical values for peak heights. The methylation status of the two genes in each sample was calculated by using the average percentages of methylation across all targeted CpG sites, respectively. For methylation analysis comparing normal adjacent and cholangiocarcinoma samples, a cutoff was calculated from tested adjacent normal samples that was based on the mean methylation value across samples + 2.069 × SD (cutoff percentage of methylation in HOXA9 = 21.4% and in HOXD9 = 18.0%). Increased methylation was defined as values above this cutoff.
Bioinformatic and statistical analysis
Unsupervised clustering.
Unsupervised hierarchical clustering was performed using DNA methylation profiles from primary cholangiocarcinoma and adjacent normal samples in TIGR MultiExperiment Viewer version 4.5 (http://mev.tm4.org) using Euclidean distance as the distance metric and average linkage clustering to merge clusters. CpG sites on the X and Y chromosomes were excluded from the analysis to avoid gender-related effects, leaving 26,486 autosomal CpG sites (14,862 CpG sites within CpG islands, 11,624 CpG sites outside of CpG islands) for hierarchical clustering.
Gene annotation.
Gene annotation was carried out using the March 2006 (NCBI/36/hg18) assembly at the UCSC database (genome.ucsc.edu/). Individual CpG sites located within a CpG island as defined by Gardiner-Garden and Frommer (25) are linked to an associated gene if the CpG island was within 2 kb distance from the respective transcription start site (26). Non-CpG island-related CpG sites are linked to a gene if they were located within 2 kb from the transcription start site.
Differential methylation analysis.
Differential methylation at CpG sites was analyzed by a linear model for microarray data (LIMMA) comparing between the group of cholangiocarcinoma samples and the group of adjacent normal samples. A false-discovery rate (FDR) of 0.05 or less was used as a cutoff for differential methylation analysis. The difference of five technical replicates of the ovarian cancer cell line PEO1 was calculated as the expected difference (null distribution). The observed difference was calculated from the grouped samples. Using bootstrap resampling method, we extracted the same number of CpG sites (n = 26,486) from the null distribution 200 times. Each time, an FDR was calculated using the formula given here. The median of FDRs was used as an estimation of the probability that CpG sites with differential methylation greater than the selected cutoff were coincidentally identified. This median FDR was set at 0.05, corresponding to ∣Δβ∣≥ 0.15:
Gene ontology analysis.
A list of genes associated with hypermethylated or hypomethylated CpG sites in cholangiocarcinoma compared with adjacent normal samples determined by LIMMA analysis with FDR of 0.05 or less and ∣Δβ∣ of 0.15 or more was used for Interpro and gene ontology (GO) analyses using the functional annotation tool in the database DAVID (http://david.abcc.ncifcrf.gov/summary.jsp). The functional terms enriched in this set of genes were examined by a modified Fisher Exact test corrected by FDR compared with the functional terms relevant to a full list of genes presented on the HumanMethylation27 BeadChip. The significance level was determined at an FDR of 0.05 or less.
Gene set enrichment analysis.
Gene set enrichment analysis (GSEA; ref. 27) was performed to determine whether 13 predefined sets of genes associated with embryonic stem cells (Supplementary Table S2) showed statistically significant, concordant methylation differences between 28 cholangiocarcinoma samples and six adjacent normal samples. For this, β values of multiple CpG sites associated with the same gene were averaged. A full set of averaged methylation data was used instead of a subset of data determined by a cutoff in GO analysis. The significance level was set for an FDR of 0.05 or less which was estimated from a 1,000 times permutation test. The GSEA package was used in this analysis.
Correlation between methods.
Correlation analysis of methylation levels of HOXA9 and HOXD9 as detected by HumanMethylation27 BeadChips and pyrosequencing was performed using Spearman correlation.
Data deposition.
Methylation profiling data are available in the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO accession number: GSE38860), and can be accessed at: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?token=zzqplqssqyeoizs&acc=GSE38860.
Results
DNA methylation profiles of primary cholangiocarcinoma and adjacent normal tissues
We analyzed the DNA methylation status of 27,578 CpG sites in 28 primary intrahepatic cholangiocarcinoma and six matched adjacent normal samples using Infinium HumanMethylation27 BeadChips. Unsupervised hierarchical clustering was performed to determine the similarity of DNA methylation profiles across primary tumors and adjacent normal tissues using CpG methylation levels (β values) of 26,486 CpG dinucleotides obtained from each sample, but excluding CpG sites on the X and Y chromosomes to avoid gender-related DNA methylation effects. Hierarchical clustering analysis returned two major clusters broadly separating between primary cholangiocarcinoma tissues and adjacent normal tissues (Fig. 1 and Supplementary Fig. S1). Cholangiocarcinoma samples formed several subclusters with one adjacent normal tissue (W022N) being embedded in the cluster of cholangiocarcinoma. Adjacent normal tissues, on the other hand, clustered tightly together with only two tumors (P055T and U030T) co-clustering with this subgroup of tissues. Reevaluation of histopathology of an adjacent biopsy to the one used to extract DNA of these exceptions confirmed normal pathology and tumor pathology, respectively. Similarly, reevaluation of tumors showed more than 80% tumor by histologic cellularity.
Unsupervised hierarchical clustering of DNA methylation profiles of primary cholangiocarcinoma (CCA) and matched adjacent normal tissues. Clustering analysis of all 26,486 autosomal CpG sites returned two major clusters, separating between primary cholangiocarcinoma tissues and adjacent normal tissues but with three exceptions: an adjacent normal sample W022N which clustered with cholangiocarcinoma, and primary cholangiocarcinoma P055T and U030T which clustered with normal.
Unsupervised hierarchical clustering of DNA methylation profiles of primary cholangiocarcinoma (CCA) and matched adjacent normal tissues. Clustering analysis of all 26,486 autosomal CpG sites returned two major clusters, separating between primary cholangiocarcinoma tissues and adjacent normal tissues but with three exceptions: an adjacent normal sample W022N which clustered with cholangiocarcinoma, and primary cholangiocarcinoma P055T and U030T which clustered with normal.
Distinctive hyper- and hypomethylation patterns in cholangiocarcinoma as compared with adjacent normal tissues
We further examined differences in methylation (∣Δβ∣ ≥ 15%) between the 28 cholangiocarcinoma and six adjacent normal tissues at all 26,486 autosomal CpG sites by using LIMMA (28). We identified 1,610 methylation changes with 809 CpG sites (associated with 603 genes) being hypermethylated and 801 CpG sites (associated with 712 genes) being hypomethylated in 28 cholangiocarcinoma compared with six adjacent normal samples at an FDR of less than 5% (Fig. 2, Supplementary Tables S3 and S4). Within the set of 28 cholangiocarcinoma, six tumors were matched to the adjacent normal samples used in this study. Among the methylation changes identified in the group comparison, 781/809 hypermethylated CpG sites (associated with 585 genes) and 748/801 hypomethylated CpG sites (associated with 665 genes) remained significant (FDR ≤ 0.05 and ∣Δβ∣ ≥ 15%) when comparing the six pairs of matched tumors and adjacent normal samples (Supplementary Tables S5 and S6).
Differential methylation in tumors (n = 28) as compared with adjacent normal samples (n = 6). The difference of methylation (x-axis) versus −log10 (FDR) (y-axis) obtained from LIMMA analysis is shown. The broken red line shows the significance cutoff corresponding to an FDR less than 0.05 and difference ∣Δβ∣ of 0.15 or more used in this study. Red spots indicate hypermethylated CpG sites; blue spots show hypomethylated CpG sites. The difference cutoff was estimated from a null distribution of differential methylation between five technical replicates on Infinium's HumanMethylation27 BeadChips.
Differential methylation in tumors (n = 28) as compared with adjacent normal samples (n = 6). The difference of methylation (x-axis) versus −log10 (FDR) (y-axis) obtained from LIMMA analysis is shown. The broken red line shows the significance cutoff corresponding to an FDR less than 0.05 and difference ∣Δβ∣ of 0.15 or more used in this study. Red spots indicate hypermethylated CpG sites; blue spots show hypomethylated CpG sites. The difference cutoff was estimated from a null distribution of differential methylation between five technical replicates on Infinium's HumanMethylation27 BeadChips.
In order to explore any patterns within the set of genes showing differential methylation at linked CpG sites in cholangiocarcinoma versus adjacent normal samples, we interrogated GO (29) and InterPro categories (30) to identify enrichment of gene functions. Table 1 shows the top 20 highly significant results for genes being hypermethylated in cholangiocarcinoma compared with adjacent normal samples with P value less than 3.62 × 10−5 after controlling for an FDR less than 0.05. The homeobox class was the most significantly enriched category within the set of genes associated with hypermethylated CpG sites in cholangiocarcinoma (FDR = 2.18 × 10−20). In addition, molecular functions and biologic processes related to homeobox, such as transcription factor activity (FDR = 1.89 × 10−11) or embryonic morphogenesis (FDR = 3.66 × 10−9) were found among the top 10 terms (Table 1). Among the set of genes which were hypomethylated in cholangiocarcinoma versus adjacent normal tissue, we observed significant enrichment of gene groups associated with the plasma membrane (FDR = 3 × 10−4), immune response (FDR = 0.005212), and the regulation and induction of apoptosis (FDR = 0.016337 and FDR = 0.032143, respectively; Table 2).
Significantly enriched top 20 functional terms in genes associated with hypermethylated CpG sites in cholangiocarcinoma
Category . | Term . | No. of genes . | % . | P valuea . | FDR . |
---|---|---|---|---|---|
INTERPRO | IPR001356:Homeobox | 46 | 8.57 | 1.18E−20 | 2.18E−20 |
INTERPRO | IPR017970:Homeobox, conserved site | 46 | 8.57 | 1.52E−20 | 2.80E−20 |
INTERPRO | IPR012287:Homeodomain-related | 46 | 8.57 | 2.16E−19 | 3.99E−19 |
GOTERM_MF | GO:0003700 ∼ transcription factor activity | 88 | 16.39 | 8.33E−12 | 1.89E−11 |
GOTERM_CC | GO:0005576 ∼ extracellular region | 126 | 23.46 | 3.07E−11 | 1.36E−10 |
GOTERM_BP | GO:0030182 ∼ neuron differentiation | 54 | 10.06 | 4.77E−10 | 3.22E−10 |
GOTERM_MF | GO:0043565 ∼ sequence-specific DNA binding | 63 | 11.73 | 2.20E−10 | 5.00E−10 |
GOTERM_BP | GO:0048598 ∼ embryonic morphogenesis | 43 | 8.01 | 5.42E−09 | 3.66E−09 |
GOTERM_BP | GO:0003002 ∼ regionalization | 34 | 6.33 | 1.25E−08 | 8.45E−09 |
GOTERM_BP | GO:0009952 ∼ anterior/posterior pattern formation | 28 | 5.21 | 3.42E−08 | 2.31E−08 |
GOTERM_BP | GO:0007389 ∼ pattern specification process | 38 | 7.08 | 9.84E−08 | 6.64E−08 |
INTERPRO | IPR001827:Homeobox protein, antennapedia type, conserved site | 13 | 2.42 | 3.27E−08 | 6.03E−08 |
GOTERM_BP | GO:0001501 ∼ skeletal system development | 41 | 7.64 | 5.28E−07 | 3.56E−07 |
GOTERM_BP | GO:0048562 ∼ embryonic organ morphogenesis | 25 | 4.66 | 8.86E−07 | 5.98E−07 |
GOTERM_BP | GO:0048568 ∼ embryonic organ development | 28 | 5.21 | 1.81E−06 | 1.22E−06 |
GOTERM_BP | GO:0048705 ∼ skeletal system morphogenesis | 23 | 4.28 | 3.86E−06 | 2.61E−06 |
GOTERM_BP | GO:0007610 ∼ behavior | 48 | 8.94 | 1.67E−05 | 1.13E−05 |
GOTERM_BP | GO:0048704 ∼ embryonic skeletal system morphogenesis | 16 | 2.98 | 2.45E−05 | 1.66E−05 |
GOTERM_BP | GO:0050877 ∼ neurologic system process | 68 | 12.66 | 3.34E−05 | 2.26E−05 |
GOTERM_BP | GO:0006355 ∼ regulation of transcription, DNA-dependent | 105 | 19.55 | 3.62E−05 | 2.45E−05 |
Category . | Term . | No. of genes . | % . | P valuea . | FDR . |
---|---|---|---|---|---|
INTERPRO | IPR001356:Homeobox | 46 | 8.57 | 1.18E−20 | 2.18E−20 |
INTERPRO | IPR017970:Homeobox, conserved site | 46 | 8.57 | 1.52E−20 | 2.80E−20 |
INTERPRO | IPR012287:Homeodomain-related | 46 | 8.57 | 2.16E−19 | 3.99E−19 |
GOTERM_MF | GO:0003700 ∼ transcription factor activity | 88 | 16.39 | 8.33E−12 | 1.89E−11 |
GOTERM_CC | GO:0005576 ∼ extracellular region | 126 | 23.46 | 3.07E−11 | 1.36E−10 |
GOTERM_BP | GO:0030182 ∼ neuron differentiation | 54 | 10.06 | 4.77E−10 | 3.22E−10 |
GOTERM_MF | GO:0043565 ∼ sequence-specific DNA binding | 63 | 11.73 | 2.20E−10 | 5.00E−10 |
GOTERM_BP | GO:0048598 ∼ embryonic morphogenesis | 43 | 8.01 | 5.42E−09 | 3.66E−09 |
GOTERM_BP | GO:0003002 ∼ regionalization | 34 | 6.33 | 1.25E−08 | 8.45E−09 |
GOTERM_BP | GO:0009952 ∼ anterior/posterior pattern formation | 28 | 5.21 | 3.42E−08 | 2.31E−08 |
GOTERM_BP | GO:0007389 ∼ pattern specification process | 38 | 7.08 | 9.84E−08 | 6.64E−08 |
INTERPRO | IPR001827:Homeobox protein, antennapedia type, conserved site | 13 | 2.42 | 3.27E−08 | 6.03E−08 |
GOTERM_BP | GO:0001501 ∼ skeletal system development | 41 | 7.64 | 5.28E−07 | 3.56E−07 |
GOTERM_BP | GO:0048562 ∼ embryonic organ morphogenesis | 25 | 4.66 | 8.86E−07 | 5.98E−07 |
GOTERM_BP | GO:0048568 ∼ embryonic organ development | 28 | 5.21 | 1.81E−06 | 1.22E−06 |
GOTERM_BP | GO:0048705 ∼ skeletal system morphogenesis | 23 | 4.28 | 3.86E−06 | 2.61E−06 |
GOTERM_BP | GO:0007610 ∼ behavior | 48 | 8.94 | 1.67E−05 | 1.13E−05 |
GOTERM_BP | GO:0048704 ∼ embryonic skeletal system morphogenesis | 16 | 2.98 | 2.45E−05 | 1.66E−05 |
GOTERM_BP | GO:0050877 ∼ neurologic system process | 68 | 12.66 | 3.34E−05 | 2.26E−05 |
GOTERM_BP | GO:0006355 ∼ regulation of transcription, DNA-dependent | 105 | 19.55 | 3.62E−05 | 2.45E−05 |
aBonferroni corrected.
Significantly enriched functional terms in genes associated with hypomethylated CpG sites in cholangiocarcinoma
Category . | Term . | No. of genes . | % . | P valuea . | FDR . |
---|---|---|---|---|---|
GOTERM_CC | GO:0044459 ∼ plasma membrane part | 123 | 20.16 | 8.27E−05 | 3.00E−04 |
GOTERM_BP | GO:0006955 ∼ immune response | 49 | 8.03 | 0.007458 | 0.005212 |
GOTERM_BP | GO:0042981 ∼ regulation of apoptosis | 53 | 8.69 | 0.023192 | 0.016337 |
GOTERM_CC | GO:0005887 ∼ integral to plasma membrane | 72 | 11.80 | 0.004246 | 0.01543 |
GOTERM_BP | GO:0043067 ∼ regulation of programmed cell death | 53 | 8.69 | 0.030145 | 0.02131 |
GOTERM_BP | GO:0010941 ∼ regulation of cell death | 53 | 8.69 | 0.033297 | 0.023575 |
GOTERM_BP | GO:0006917 ∼ induction of apoptosis | 28 | 4.59 | 0.045123 | 0.032143 |
GOTERM_BP | GO:0012502 ∼ induction of programmed cell death | 28 | 4.59 | 0.047725 | 0.034042 |
GOTERM_CC | GO:0031226 ∼ intrinsic to plasma membrane | 72 | 11.80 | 0.009136 | 0.033281 |
Category . | Term . | No. of genes . | % . | P valuea . | FDR . |
---|---|---|---|---|---|
GOTERM_CC | GO:0044459 ∼ plasma membrane part | 123 | 20.16 | 8.27E−05 | 3.00E−04 |
GOTERM_BP | GO:0006955 ∼ immune response | 49 | 8.03 | 0.007458 | 0.005212 |
GOTERM_BP | GO:0042981 ∼ regulation of apoptosis | 53 | 8.69 | 0.023192 | 0.016337 |
GOTERM_CC | GO:0005887 ∼ integral to plasma membrane | 72 | 11.80 | 0.004246 | 0.01543 |
GOTERM_BP | GO:0043067 ∼ regulation of programmed cell death | 53 | 8.69 | 0.030145 | 0.02131 |
GOTERM_BP | GO:0010941 ∼ regulation of cell death | 53 | 8.69 | 0.033297 | 0.023575 |
GOTERM_BP | GO:0006917 ∼ induction of apoptosis | 28 | 4.59 | 0.045123 | 0.032143 |
GOTERM_BP | GO:0012502 ∼ induction of programmed cell death | 28 | 4.59 | 0.047725 | 0.034042 |
GOTERM_CC | GO:0031226 ∼ intrinsic to plasma membrane | 72 | 11.80 | 0.009136 | 0.033281 |
aBonferroni corrected.
We validated our array-based findings in a bigger panel of cholangiocarcinoma (n = 102) and adjacent normal (n = 24) samples using bisulfite pyrosequencing to analyze differentially methylated loci at the two HOX genes—HOXA9 and HOXD9—identified as differentially methylated. Hypermethylation of HOXA9 and HOXD9 was, respectively, found in 86.3% (88 out of 102) and 89.2% (91 out of 102) of primary cholangiocarcinoma samples. Primary cholangiocarcinoma showed significantly increased methylation (P < 0.001, unpaired t test) as compared with adjacent normal samples at both loci (Fig. 3A and B). We found a strong correlation between pyrosequencing-based detection of CpG methylation and β values obtained from HumanMethylation27 BeadChips (Spearman correlation coefficient > 0.9, P < 0.001) for both genes supporting the observed hypermethylation of linked CpG sites in cholangiocarcinoma.
Differential methylation at loci on HOXA9 (A) and HOXD9 (B) in primary cholangiocarcinoma (CCA) compared with adjacent normal samples as detected by bisulfite pyrosequencing. Scatter-dot plots represent percentage of methylation, and black horizontal lines are average methylation levels across adjacent normal and primary cholangiocarcinoma samples, respectively. Significant differences of methylation levels are observed between primary cholangiocarcinoma, compared with adjacent normal samples (P < 0.001, unpaired t test, GraphPad Prism software).
Differential methylation at loci on HOXA9 (A) and HOXD9 (B) in primary cholangiocarcinoma (CCA) compared with adjacent normal samples as detected by bisulfite pyrosequencing. Scatter-dot plots represent percentage of methylation, and black horizontal lines are average methylation levels across adjacent normal and primary cholangiocarcinoma samples, respectively. Significant differences of methylation levels are observed between primary cholangiocarcinoma, compared with adjacent normal samples (P < 0.001, unpaired t test, GraphPad Prism software).
Hypermethylation of PRC2 targets is enriched in cholangiocarcinoma compared with adjacent normal tissue
Homeobox genes have a vital role in embryonic development and differentiation. Regulation of HOX involves the epigenetic key mark trimethylation at histone H3 lysine 27 (H3K27me3), which is established and maintained by the polycomb group of proteins (PcG), particularly the polycomb repressive complex 2 (PRC2) which maintains early patterns of HOX gene repression (31). It has been previously suggested that PRC2 target genes in embryonic stem cells may be more vulnerable to epigenetic silencing by DNA methylation. PcG protein repression targets many more genes in addition to the homeobox genes (32, 33) in undifferentiated embryonic cells; therefore, we addressed whether additional targets of PcG proteins would be enriched within the set of aberrantly methylated genes in cholangiocarcinoma. To this end, we performed GSEA interrogating 13 published embryonic gene signature sets associated with human embryonic stem cells (hES) for their methylation pattern at loci associated with CpG islands in cholangiocarcinoma and adjacent normal tissues (Supplementary Table S2). Five signatures were significantly hypermethylated in cholangiocarcinoma as compared with adjacent normal at P values less than 0.05 and an FDR less than 0.05 encompassing target genes of the PRC2 components EED (P = 0.009 and FDR = 0.04) and SUZ12 (P = 0.017 and FDR = 0.04), targets of H3K27me3 (P = 0.017 and FDR = 0.04), PRC2 target genes (shared targets of EED, SUZ12, and H3K27me3; P = 0.011 and FDR = 0.04) and homeobox genes (P = 0.014 and FDR = 0.04; Table 3). Overall, the observed enrichment of differential methylation in cholangiocarcinoma at embryonic gene sets, particularly those associated with the PRC2 and its components, indicates aberrant DNA methylation of genes associated with a more stem cell-like phenotype in cholangiocarcinoma tumors as compared with adjacent normal tissue.
GSEA between cholangiocarcinoma (n = 28) and adjacent normal (n = 6) samples
Methylation status . | Gene seta . | GSA Score . | P value . | FDR . |
---|---|---|---|---|
Hypermethylation | EED targets | 0.2396 | 0.009 | 0.04b |
PRC2 targets (1) | 0.3557 | 0.011 | 0.04b | |
H3K27me3 bound | 0.211 | 0.017 | 0.04b | |
Homeobox genes | 0.6098 | 0.014 | 0.04b | |
SUZ12 targets | 0.2559 | 0.017 | 0.04b | |
NOSTF targets | 0.1956 | 0.158 | 0.34 | |
Hypomethylation | PRC2 targets (2) | −0.402 | 0.045 | 0.33 |
OCT4 targets | −0.1652 | 0.093 | 0.33 | |
NOS targets | −0.1803 | 0.106 | 0.33 | |
NANOG targets | −0.1923 | 0.125 | 0.33 | |
ES2 | −0.2615 | 0.125 | 0.33 |
Methylation status . | Gene seta . | GSA Score . | P value . | FDR . |
---|---|---|---|---|
Hypermethylation | EED targets | 0.2396 | 0.009 | 0.04b |
PRC2 targets (1) | 0.3557 | 0.011 | 0.04b | |
H3K27me3 bound | 0.211 | 0.017 | 0.04b | |
Homeobox genes | 0.6098 | 0.014 | 0.04b | |
SUZ12 targets | 0.2559 | 0.017 | 0.04b | |
NOSTF targets | 0.1956 | 0.158 | 0.34 | |
Hypomethylation | PRC2 targets (2) | −0.402 | 0.045 | 0.33 |
OCT4 targets | −0.1652 | 0.093 | 0.33 | |
NOS targets | −0.1803 | 0.106 | 0.33 | |
NANOG targets | −0.1923 | 0.125 | 0.33 | |
ES2 | −0.2615 | 0.125 | 0.33 |
aOnly gene sets with FDR < 50% are shown in the table.
bFDR ≤ 0.05.
Discussion
Aberrant DNA methylomes are one of the key features of cancer (11). Recent advances in DNA methylation technology have allowed for genome-wide profiling of cancer tissues providing a vast amount of previously unknown methylation changes that furthers our understanding of the underlying pathways associated with tumorigenesis. To our knowledge, this is the first report on a genome-wide DNA methylome analysis in cholangiocarcinoma, a cancer associated with liver-fluke infection in endemic countries. We show that DNA methylation is a common molecular aberration in cholangiocarcinoma, with DNA methylation profiles of primary cholangiocarcinoma being different from those of adjacent normal tissues. Differential methylation between cholangiocarcinoma and adjacent normal tissues comprised hypermethylation and hypomethylation events which were equally present at CpG sites within and outside of CpG islands. As previously discussed, our findings of distinct methylation changes in cholangiocarcinoma as compared with normal tissues relate to liver-fluke–associated cholangiocarcinoma, whereas the sporadic type found in Western countries has not been included in this study. Although cholangiocarcinoma from different geographic sources may share common biologic mechanisms occurring during tumorigenesis, additional studies are required to analyze differential methylation in Western patients with cholangiocarcinoma versus normal tissues and to filter out differences between Western and Thai patients with cholangiocarcinoma. To demonstrate the utility, in relevant biologic processes, of differential methylation between cholangiocarcinoma and adjacent normal samples, we further showed that hypermethylation in cholangiocarcinoma is significantly overrepresented at genes involved in the developmental processes, most importantly homeobox and genes involved in transcription factor activity and embryonic morphogenesis. More specifically, we examined methylation states of gene sets associated with a stem-cell like phenotype, and identified hypermethylation at specific target genes of the PcG complexes including PRC2 targets, homeobox genes, EED targets, SUZ12 targets, and H3K27me3 target genes. In addition, more than 55% of homeobox genes are the target genes of PRC2, EED, SUZ12, and H3K27me3, indicating that the overlap of these gene sets might be the main target of hypermethylation in cholangiocarcinoma. PcG complexes are involved in transcriptional repression of genes. They play a role to maintain heritable transcription patterns of the homeobox genes during development and differentiation, in conjunction with their positively acting counterparts known as the Trithorax group (TrxG) of proteins (34, 35). A previous report has described that stem cell PcG targets are up to 12-fold more likely to have cancer-specific promoter DNA hypermethylation than nontargets. These genes might be permanently repressed leading to loss of differentiation properties of the cells (36). Several reports have described that tumors often lack the terminal differentiation traits possessed by their normal cell states. These parallels have supported the hypothesis that cancer cells arise from undifferentiated stem or progenitor cells or, alternatively, that cancer cells can undergo progressive dedifferentiation during their development (21–23). Taken together, our findings support the theory of a stem cell origin for cholangiocarcinoma.
Consistent with the observation of a stem-cell associated phenotype and in support of our observation of hypermethylation of PRC2 target genes in cholangiocarcinoma, a recent study has shown overexpression of EZH2, the catalytic subunit of PRC2, in cholangiocarcinoma (37). Aberrant levels of EZH2 have previously been linked to tumor development and have been observed in a variety of tumors including prostate, breast, and melanoma (38–40). Moreover, high EZH2 levels have been associated with poor prognosis in these tumors types. In cholangiocarcinoma, EZH2 expression levels were significantly and positively associated with step-wise progression of cholangiocarcinoma from low-grade dysplasia, high-grade dysplasia, carcinoma in situ, to invasive carcinoma, suggesting that chromatin remodeling complexes are involved in the development and progression of this cancer type (37). High EZH2 was linked to increased methylation of the p16Ink4a promoter (37), supporting the notion that alternative complementary silencing mechanisms such as DNA methylation may induce or permanently replace initial silencing of PRC2 target genes. In our study, low methylation levels of the homeobox genes HOXA9 and HOXD9, both PRC2 target genes, were verified in adjacent normal tissues. Although the number of adjacent normal tissues initially used for DNA methylome analysis was comparatively small, the inclusion of more normal adjacent samples in our pyrosequencing analyses indicates that low methylation levels can be observed across additional normal samples at homeobox genes. Methylation of homeobox in normal tissue could reflect a mark of rare cancer stem cells, which then spread in tumors. In this respect, our observation of DNA hypermethylation at PRC2 target genes including homeobox in cholangiocarcinoma shows potential to serve as a biomarker for early detection of the disease. In addition, PRC2 and EZH2, in particular, might further be attractive therapeutic and, even more important, chemopreventative targets in cholangiocarcinoma. In glioblastoma multiforme, an aggressive type of brain tumor, pharmacologic EZH2 inhibitors such as DZNep have been shown to strongly impair cancer stem cell self-renewal in vitro and tumor-initiating capacity in vivo (41). Therefore, DZNep as well as novel, more specific pharmacologic agents (42) targeting histone lysine methyltransferases may have potential to serve as a novel approach for early intervention of the often terminal cholangiocarcinoma.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Authors' Contributions
Conception and design: R. Sriraksa, R. Brown
Development of methodology: R. Sriraksa, W. Dai, R. Brown
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): C. Zeller, A. Siddiq, A.J. Walley, T. Limpaiboon, R. Brown
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): R. Sriraksa, C. Zeller, W. Dai, R. Brown
Writing, review, and/or revision of the manuscript: R. Sriraksa, C. Zeller, A.J. Walley, T. Limpaiboon, R. Brown
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): R. Sriraksa
Study supervision: T. Limpaiboon, R. Brown
Acknowledgments
The authors thank the Liver Fluke and Cholangiocarcinoma Research Center, Khon Kaen University, for providing samples and clinical data.
Grant Support
This work was financially supported by the Higher Education Research Promotion and National Research University Project of Thailand, the Office of the Higher Education Commission, through the Center of Excellence in Specific Health Problems in the Greater Mekong Sub-region cluster (SHeP-GMS), Khon Kaen University, Thailand (grant nos. PD54201 and NRU542015 to R. Sriraksa and T. Limpaiboon, respectively), and by the Imperial Experimental Cancer Medicine Centre and Cancer Research UK program (grant no. C536/A13086 to R. Brown).
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.