Cancer associated epigenetic transitions identified by genome-wide histone methylation binding profiles in human colorectal cancer samples and paired normal mucosa

H3K4me3 and H3K27me3 binding profiles in paired normal and tumor samples from CRC patients

We generated enrichment profiles for H3K4me3 and H3K27me3 in tumor samples from two CRC patients, for which the profiles in normal mucosa have previously been generated [13] and for the colon adenocarcinoma cell line HT29. After determining the enriched regions (Figure 1A) for each modification in normal and tumor samples (Table 1 and Additional file 1), they were annotated to human genes using the UCSC knownGenes databases [19] (Table 1) (Methods).

Signal footprints were created around transcription start sites and both the H3K4me3 and H3K27me3 patterns were similar to those reported [20–22] in other cell types (Figure 1B). As a quality control, microarray expression data from public repositories for normal colon, CRC tumors and HT29 were used [23]. Although the expression data comes from unrelated patient material, the objective was to study genome wide trends in transcription given the epigenetic status and not single gene effects where expression from the same patient material is essential. The expression of H3K4me3 enriched genes were higher than average, while H3K27me3 genes tended to be silent (Figure 2A). Bivalent and especially semi-bivalent (Methods) gene groups included fewer genes and were poorly expressed, although with higher variation in expression. The quality of our data is also supported by previous validation of H3K4me3 and H3K27me3 enriched regions by qPCR in normal colon [13] demonstrating that ChIP-chip analysis can be performed using human tissue as biological material.

For regions with associated transcripts, gene ontology (GO) annotations were obtained and an analysis was performed for over/under representation of terms compared to a whole genome background (Figure 2B, Methods). The H3K4me3 and H3K27me3 gene groups have opposite functions and bivalent and semi-bivalent gene groups lie closer in functional space to H3K27me3 than H3K4me3 genes, as also suggested by expression data. Therefore, bivalent and semi-bivalent genes were treated as part of H3K27me3 gene groups for most subsequent analysis. In general, H3K4me3 fractions displayed over-representation in metabolic processes while H3K27me3 fractions showed enrichment in developmental categories. The overrepresented categories in H3K4me3 and H3K27me3 genes were very similar in the HT29 cell line.

H3K4me3 and H3K27me3 profiles are similar between normal and tumor samples but distinct compared to CRC cell lines

Using all probes we calculated the Pearson's correlation between all combinations of histone modifications, patient material (normal/tumor) and cell line (Additional file 2, Table S1). In all cases the correlation between any patient materials was always strongest when that material was compared to another patient material. In some cases, e.g. patient 1 tumor H3K27me3/K4me3, the correlation was highest with patient 2 normal H3K27me3/K4me3 rather than the HT29 cancer cell line. This pattern is more evident when regarding enriched regions only and by counting overlapping enriched regions, a relatively high similarity between normal and tumor samples and between patients for both H3K4me3 (Figure 3) and H3K27me3 was noted (Table 1 and Additional file 2, Figure S6). Around 60% of the H3K4me3 and 35% of H3K27me3 regions in both patients are common between normal and tumor tissue. When annotated transcripts rather than individual regions were used these numbers rise to around 70% and 50% respectively. Considering H3K4me3 genes, HT29 targets showed large overlaps with both patients normal and tumor samples (Table 1-2). The overlaps for the adenocarcinoma cell line SW480 [24] were quite poor (10-20%) both with tissue samples and HT29 (Table 2). We also noted that the number (around 1%, 80-431) of bivalent regions detected in promoters here is lower than the number detected in a previous study [25] on differentiated cells (CDC36+) where around 3% (693) of the genes were found to have bivalently marked promoters. Cui et al [25] performed their study using next generation sequencing and called more regions enriched throughout the genome suggesting that the main cause of difference lies in stringency of the downstream analysis.

Detection of genes and pathways with chromatin state differences between tumor and normal samples

To identify the genomic regions differing in their chromatin state between normal and tumor tissues (Methods), the results were split in fifteen so called transitions each representing a combined state of one or more marks in normal and tumor. Three major groups were created; i) gaining or loosing a single mark, ii) keeping a mark and iii) changes involving both marks (Figure 4A and Additional file 2, Figure S4A). The majority fall into category ii), that is no change, and category i) which is in agreement with the large overlaps reported in the previous section (Additional file 2,Table S2).

Category iii) implies chromatin transitions between normal and tumor that can result in a switch from transcriptionally active (e.g. H3K4me3) to inactive states (e.g. H3K27me3) or vice verse. Therefore, it is the most interesting category in terms of potentially functional consequences from an oncogenic point of view. However, only 4 genes were common to both patients in category iii) (Additional file 2, Table S2). KLF7 (N K27 to T K4), a transcription factor belonging to the KLF family that plays critical roles in differentiation, development, and maintenance of tissue homeostasis [26], but, to the best of our knowledge, without a specific role in colon cancer. The remaining three common genes (EBF3, DKFZp667I0324 and RBMS1) were all bivalent in normal tissue but lost all their histone methylation in tumor. Among these genes, 2 have previously been functionally characterized. EBF3 is a transcription factor reported to be a tumor suppressor gene [27] with observed silencing in CRC cell line HCT116, was bivalent in normal and lost both marks in tumor. RBMS1 is a transcription factor that has been suggested as an important factor in inducing apoptosis in mice [28], and has been shown to have predictive power on colorectal cancer recurrence in human patient cohorts [29]. A number of transitions were specific to individual patients e.g. from category iii). Some genes had H3K4me3 in normal colon which changed to H3K27me3 in tumors, e.g. in patient 1 HOXB13 and KREMEN2, which are frequently down regulated in CRC [30, 31]. Among genes with a silent state (H3K27me3 mark) in normal colon that became active in CRC samples (gained H3K4me3) there were several well-characterized oncogenes [32–36], e.g. in patient 1 FLI1, WWTR1 & ZEB2 and in patient 2 TACSTD2 & TWIST1 (Figure 1A and Additional file 2, Figure S5 (patient 1)), the last two with proposed importance in CRC [32, 35].

Several genes in categories i) and ii) were found to have the same transitions in both patients. Both our patients have H3K4me3 on the PTGER2 gene in normal tissue and gained H3K27me3 in the tumor, which offers a mechanistic explanation to the previously reported expression pattern [37] of this gene in colorectal cancer. In addition, three genes common to both patients which carried bivalent marks in normal and lost H3K4me3 in tumor, LHX9, PKNOX2 and LBXCOR1 and one, PRDM8, lost the H3K27me3 in tumor, all of which are transcription factors with unknown function in colon.

The chromatin states of genes involved in biological pathways with relevance in CRC biology [38] was specifically investigated. The WNT-signalling pathway has been implicated in CRC [39, 40], with the APC gene frequently silenced or mutated. Patient 1 has tumor repressive transitions in the APC gene and tumor activating transitions in Axin, PP2A, β-catenin and in many downstream factors along the WNT-signalling pathway affecting transcription, cell cycle and CRC development, especially within the TCF/LEF family [41, 42] (Figure 4B). This suggests a repression on the β-catenin degradation complex and an up regulation of β-catenin dependent gene expression with downstream effects on the cell cycle [43]. Recently, a dual role of APC in the WNT-signalling pathway has been suggested with both activating and repressing mechanisms [44] which could explain the patterns found in patient 2 where APC has gained activating markers (Additional file 2, Figure S4B) In addition, patient 1 shows an activation of several known repressors of the WNT-pathway. CtBP represses E-cadherin which in turn is a tumor repressor restricting tumor cell motility and invasion [45]. The ubiquitin ligase Cul1 is part of a mechanism targeting proteins for degradation and thus an important factor in cell cycle control. Specifically, Cul1 acts as an oncogene that forms complexes that specifically target the cell cycle inhibitor p27/Kip1 for degradation [46].

Global expression changes of H3K4me3 and H3K27me3 targets in CRC

Although epigenetic transitions within category iii) are expected to have important functional consequences, few examples are detected in a given CRC patient. Therefore, these examples are less likely to have large-scale changes in gene expression. More global relationships between chromatin state and gene expression were analyzed using all epigenetic transitions involving more than 100 genes i.e. maintenance of H3K4me3 or H3K27me3, tumor loss of H3K4me3 or H3K27me3 as well as tumor gain of H3K4me3 or H3K27me3. Microarray expression data for paired normal and tumor samples from 24 CRC patients [23] were used, and the expression changes between normal and cancer samples for the genes included in each transition category (patient specific) was compared (Figure 5 and Additional file 2, Figures S1-S3).

Genes in transitions where H3K4me3 was maintained or specifically gained in tumors, showed significantly higher expression levels in cancer samples (Figure 5A and Additional file 2, Figure S2). Opposite expression shifts were observed for genes that kept, gained or lost H3K27me3 in the tumor samples (Figure 5B, Additional file 2, Figure S2). The H3K4me3 and H3K27me3 gene groups also differed in the baseline expression in normal samples as H3K4me3 genes were highly expressed, while most H3K27me3 genes showed low expression levels. Furthermore, genes in H3K4me3 transition groups are typically already activated in normal tissue and become "hyperactivated" in tumors, while genes in H3K27me3 transition groups are lowly expressed in normal colon and become "hypersilenced" in cancer samples. The generality of the observations is supported by (i) the similar expression shifts obtained when considering the epigenetic transition groups from either patient 1 and patient 2 (Additional file 2, Figure S2A); and (ii) from a second expression data set with paired normal and tumor samples from 9 CRC obtained in a different array platform [47] (Additional file 2, Figure S2 B-C).

The significantly up- or down- regulated genes in cancer samples in the two expression microarray data sets described above assuming a paired sample design was analyzed further. For both patients and microarray data sets (Figure 5C and Additional file 2, Figure S3) genes keeping or gaining H3K4me3 in tumors were significantly over-represented among up-regulated genes and significantly under-represented among down-regulated genes. Opposite patterns were observed for genes that kept, gained or lost H3K27me3 in tumor samples (Figure 5C and, Additional file 2, Figure S3). These results indicate that at least for a fraction of genes the reported expression changes may be of large magnitude and occur in most patients. We gathered microarray expression data for HT29 [23] and considered the same epigenetic transitions as described above, and can report that the tendencies towards hyperactivation or hypersilencing of H3K4me3 and H3K27me3 gene groups became even more pronounced in HT29 cells (Figure 5D-E and Additional file 2, Figure S4).

Hypersilencing of some H3K27me3 genes may involve loss of H3K27me3 and DNA hypermethylation

The expression levels between normal colon, tumor samples and HT29 cell line for all H3K4me3 and H3K27me3 enriched genes were compared in the three samples types and for those genes that were shared or unique to any of them. The overall H3K4me3 or H3K27me3 gene groups showed progressive hyperactivation and hypersilencing in tumor and HT29 cells compared to normal colon, respectively (Figure 6A,B). This was also true for the genes with shared histone modification states, but less obvious in some of the unique groups. For the H3K27me3 unique groups, the hypersilencing pattern was still found for those genes carrying that histone mark only in normal colon samples (Figure 6B).

It has been recently reported [12] that genes enriched in H3K27me3 in normal prostate frequently become DNA hypermethylated and lose H3K27me3 in prostate cancer cell lines. Furthermore, DNA hypermethylation seems to be exacerbated by in vitro culture [14]. In order to test whether genes with the H3K27me3 mark in normal colon overlapped more with genes hypermethylated in CRC than genes that carried the histone mark in tumor samples, the lists of H3K27me3 genes was compared with those previously reported as frequently hypermethylated in CRC [48]. As shown in Figure 6C, both H3K27me3 genes in normal and tumor colon samples were significantly over-represented among DNA hypermethylated genes. This was not the case for the HT29 genes. Over-representation was only significant for normal colon, when considering H3K27me3 unique groups, and in fact even a trend towards under-representation was observed for H3K27me3 genes unique to HT29 (Figure 6C).

Patient material

Normal colon mucosa [13] and tumor samples were obtained from two patients diagnosed with CRC at Uppsala University Hospital (Sweden). Immediately upon surgery, the removed material was kept on ice. Samples from the normal colon mucosa were taken at least 15 cm from the localized tumor, while the tumor sample was collected without containing parts of the normal mucosa. A gastrointestinal pathologist performed sample collection and tumor evaluation. The patients did not receive any anticancer treatment before surgery. All patients treated for CRC in Uppsala are registered in the Swedish Colorectal Cancer Registry. When admitted to the hospital they are asked to participate in that registry and also informed that biopsies will be saved in a biobank. After having read the information the informed consent is verbally obtained. Linked to this we have an ethical approval to perform studies like this, #2006/077.

ChIP and ChIP-chip

ChIP on patient samples and HT29 cells were performed as previously described [13, 54]. DNA amplification, fragmentation, labelling and hybridizations of ChIP and input DNAs were performed according to Affymetrix recommendations and basically as previously described [13], using Affymetrix GeneChip Human Promoter 1.0 arrays, which cover approximately from 7.5 kb upstream to 2.45 kb downstream of transcription start sites for over 25,500 human promoters. Raw array data has been deposited in ArrayExpress under accession number E-TABM-533

Data analysis pipeline

Affymetrix Tiling Analysis SDK revision 4 was used to quantile normalize each pair of replicate ChIP measurements together with the two input measurements and adjusted to have a median intensity of 200 [55, 56]. Replicates were combined using a sliding window of 301 bp (bandwidth = 150) assigning the median value of the replicates in the window to the center probe. From this point, all analysis was carried out using custom scripts in R [57]. Log2-ratios of merged ChIP measurements over merged input measurements were used for further analysis. Enriched regions were defined using a Z-score based sliding window approach. For each probe on the array, a Z-score was calculated from the average of log2-ratios of probes within a window of 150 base pairs centred on the probe against the distribution of log2-ratios on the array.

Enriched regions were defined from enriched probes (Z-score >= 6) using the following criteria; (i) at least 2 enriched probes within the region, (ii) each enriched probe must be flanked on both sides by existing measurements within 110 base pairs (iii) the maximum gap between enriched probes within a region is 110 base pairs. Regions were prolonged with probes positioned within 75 base pairs of the probes defining the tails of the regions. Finally, regions were merged if the shortest distance was less than 1000 base pairs.

Bivalent regions were defined as regions sharing at least three probes between H3K4me3 and H3K27me3 regions out of which at least two were enriched. Regions classified as having changed their methylation status ('transitions') between normal and tumor where required to have at least three probes in common out of which two needed to be commonly labelled as enriched in the different states. In addition to this we required a Z-score of at least 3 between the normal and tumor in the non-changed methylation status ensuring that we are not selecting regions immediately on opposite sides of the chosen cut-off, i.e. clearly separating changes involving bivalent marks to single marks from instances where a single mark is changed to another single mark.

Gene Annotation

We associated enriched regions to UCSC knownGene transcripts if; (i) the regions were intragenic and positioned within 2.5 kb downstream of TSS or (ii) intergenic and positioned within 2.5 kb upstream of TSS but not associated to any other transcript through intragenic positioning. This window size (+/- 2.5 kb) was chosen based on the distribution of both H3K4me3 and H3K27me3 signals around transcription start sites. We defined transcripts to be marked by H3K4me3 if associated with at least one H3K4me3 region and no H3K27me3 region. The same was done for H3K27me3 associated transcripts. Transcripts associated with both H3K4me3 regions and H3K27me3 regions but not with bivalent regions were called semi-bivalent.

Gene Ontology (GO) analysis

GO terms were tested using a one-sided Fishers exact test. P-values were corrected for multiple hypotheses testing using Bonferroni correction. The hierarchical clustering was done using Euclidean distances and the Ward agglomeration method. All analysis was done using R [57]

Microarray Gene expression data sets

Microarray genes expression data sets were downloaded from Gene Expression Omnibus (GEO) public repositories. Expression data was obtained and analyzed using GEOquery and the 'limma' Bioconductor package in R [57]. Expression data from paired tumor and normal colon mucosa samples from 24 and 9 different CRC patients corresponded to GSE10950 and GSE5364 datasets, respectively. HT29 expression data (GSM277543) was originally generated in parallel with the tissue samples from GSE10950 dataset.