Genome-wide profiling in colorectal cancer identifies PHF19 and TBC1D16 as oncogenic super enhancers

Colorectal cancer is one of the most common cancers in the world. Although genomic mutations and single nucleotide polymorphisms have been extensively studied, the epigenomic status in colorectal cancer patient tissues remains elusive. Here, together with genomic and transcriptomic analysis, we use ChIP-Seq to profile active enhancers at the genome wide level in colorectal cancer paired patient tissues (tumor and adjacent tissues from the same patients). In total, we sequence 73 pairs of colorectal cancer tissues and generate 147 H3K27ac ChIP-Seq, 144 RNA-Seq, 147 whole genome sequencing and 86 H3K4me3 ChIP-Seq samples. Our analysis identifies 5590 gain and 1100 lost variant enhancer loci in colorectal cancer, and 334 gain and 121 lost variant super enhancer loci. Multiple key transcription factors in colorectal cancer are predicted with motif analysis and core regulatory circuitry analysis. Further experiments verify the function of the super enhancers governing PHF19 and TBC1D16 in regulating colorectal cancer tumorigenesis, and KLF3 is identified as an oncogenic transcription factor in colorectal cancer. Taken together, our work provides an important epigenomic resource and functional factors for epigenetic studies in colorectal cancer.

The binding of transcription factors (TFs) to enhancers is one of the critical steps in transcription activation. Recently, the development of epigenomics revealed novel features of active and silent enhancers and shed light on the study of transcriptional regulation in multiple research fields1,2,3,4,5. Epigenetic marks on chromatin are important signatures for cell identification, which co-operate with transcription factors to regulate transcription2,6,7. Histone modifications mark enhancers on chromatin and are critical for their activity. H3K4me1 is the mark for primed enhancers;7,8 H3K27ac for active enhancers and H3K27me3 for poised enhancers1. Though the initial discovery was concluded from ChIP-Seq of mediator subunits, now H3K27ac in the intergenic chromatin is widely used for the identification of active enhancers6,9,10. Moreover, it was discovered that many genes are often regulated by multiple enhancers and the state of these enhancers varies in different cell types1,4. Therefore, it has emerged as critical questions for many fields how enhancer activity is regulated for signaling pathways and selective gene transcription.

Pioneer studies hypothesized that gain of enhancer activity is one of the common features for cancers6,11,12,13, which is supported by some recent studies in patients and animal models10,14,15. However, it is still not clear whether it is a common feature for all cancers or just a portion of them. Interestingly, many genes related with epigenetic regulation of enhancer activity are frequently mutated in cancer, such as lysine methyltransferase 2 C/D(KMT2C/D, also named as MLL3/4), E1A binding protein p300 (EP300), CREB binding protein(CEBBP), lysine demethylase 6 A (KDM6A, also named as UTX) and lysine demethylase 5 C(KDM5C)12,16,17,18,19,20. Moreover, inhibitors for bromodomain-containing 4 (BRD4), one reader of H3K27ac on enhancer, were shown to be effective in cancer treatment21. It is then urgent to clarify the roles of enhancers in cancer and the underlying mechanisms.

It has been shown that the enhancers controlling the transcription of key oncogenic genes, such as MYC proto-oncogene (MYC), distinguish in different types of cancers22. It is probably that transcription factors activated by signaling networks vary a lot in different cancer cells, which causes that enhancers are activated in different ways in response to variant genome mutations and upstream signals. Thus, enhancer profiling may reflect the feature of distinguished cancers and be used for classification10,13,23,24,25. Genome-wide profiling of active enhancers has been carried out in several types of cancers, such as pancreatic cancer, nasopharyngeal carcinoma, and clear cell renal carcinoma, usually using the approach of H3K27ac ChIP-Seq10,14,26. A pan-cancer study using TCGA RNA-Seq data also analyzed global enhancer distribution24, but since a large portion of enhancer RNA is difficult to detect, the approach is not very reliable.

Colorectal cancer (CRC) is one of the most common cancers in the world. Recent studies about aberrant DNA methylation have gained sight of the diagnosis community, and epigenetic regulation becomes one of the critical regulatory factors for CRC27,28,29,30,31,32. Some groups have studied the genome-wide distribution of active enhancers in CRC23,33. The early studies used H3K4me1 as a mark which was not suitable to identify functional active enhancers23. H3K4me1 marks primed enhancers, which does not represent functional enhancers in the genome. A recent study took normal colonic epithelial crypts, CRC cell lines, and four primary patient colorectal tumors, and analyzed the genome-wide difference of active enhancers using H3K27ac ChIP-Seq33. Flebbe et al. also performed H3K27ac in a small number of rectal cancer tissues, but did not analyzed cancer-specific enhancer features34. These studies used very few clinical samples, which could not reflect the real clinical features of CRC tissues and are not very helpful to enhancer studies in CRC.

In this work, to establish a comprehensive map for active enhancers in CRC, we perform H3K27ac ChIP-seq analysis with 73 pairs of CRC tissues (tumor tissues with paired adjacent native tissues), as well as the corresponding genomic and transcriptomic sequencing. We identify thousands of enhancers and multiple TFs involved in CRC, and a portion of them are experimentally verified, which provides important epigenomic resources and research candidates for future studies in CRC.

Genome-wide study of enhancer distribution in CRC patient tissues

To establish a comprehensive genome-wide view of active enhancers of CRC patient tissues, we totally collected 80 pairs of tissues (tumor and their adjacent tissues) from CRC patients. The tissues were collected from patients who received surgical treatment at Zhongnan Hospital of Wuhan University (Wuhan, China), and no specific criteria were applied when collecting tissues. The patients were mostly from the Huazhong area of China, especially Hubei province. We optimized the ChIP-Seq protocol and performed H3K27ac ChIP-Seq for these samples, as well as the corresponding mRNA and input DNA sequencing. Some samples failed in the study, and eventually, we got high quality of sequencing data from 74 CRC tissues and 73 native tissues, among which 73 were paired (Fig. 1A, Supplementary Fig. S1 & Supplementary Data 1). We also performed H3K4me3 ChIP-Seq for 43 pairs of tissues. Totally we generated 524 high-quality sequencing samples, including 147 H3K27ac, 86 H3K4me3, 144 RNA-Seq, and 147 genomic sequencing samples (Sup. Data 2), which provide important epigenomic information for CRC studies.

A Experimental workflow for studying the enhancer landscapes of tumor and native tissues from CRC patients. B Genomic distribution of enhancer elements in tumor and native tissues from CRC patients. C Saturation analysis showing the percentage of newly gained enhancers comparing with total significant enhancers along with an increasing number of the tumor samples. D Fold change (FC) and p.adj of human gene expression comparing tumor and native tissues. Red dots represent tumor up-regulated genes, blue dots for native tissue up-regulated genes, and grey dots for genes not changed. E Normalized ChIP-seq and RNA-seq Meta tracks showing H3K27ac and mRNA signal on MYC promoter and enhancer loci. F Overlap of enhancer loci between our patient data and 20 COAD cell lines (GSE77737, Andrea J. Cohen et al.). G Percentage of novel enhancers in CRC identified in our study.

Full size image

For bioinformatic analysis, we first identified active enhancers in each tissue by H3K27ac peaks far away from transcription start sites (TSS). To avoid interference by artificial results, the enhancers at least appearing in two samples were considered as significant enhancers. Our analysis revealed totally of 27,156 significant enhancers in native tissues and 39,207 in tumor tissues, most of which were distributed on introns and intergenic regions as expected (Fig. 1B & Supplementary Fig. S2A). Meanwhile, we identified 9896 and 10663 active promoters in native and tumor tissues, respectively (Supplementary Fig. S2A). The saturation analysis showed that the gained enhancers in tumors reached 80% when using less than 40 pairs of samples for analysis, and 90% with around 50 pairs, indicating the sample size used in our study was good enough for statistical analysis (Fig. 1C). We downloaded H3K4me1 (ENCODE, ENCFF557VIT) and BRD4 (GEO, GSM3593876) ChIP-seq data in HCT116, a CRC cell line, from public databases, and calculated the RPM values of H3K4me1 and BRD4 signal in significant CRC tumor enhancer loci. Our analysis showed that the H3K27ac peaks of our study are nicely correlated with the published BRD4 and H3K4me1 signal (Sup. Fig. S2B&C). Our RNA-Seq analysis identified 2226 up-regulated differently expressed genes (DEGs) and 1979 down-regulated DEGs in CRC tumors (Fig. 1D). Compared with the TCGA data, many DEGs are overlapped (Sup. Fig. S2D). The proportion of genes assigned with multiple enhancers was shown (Sup. Fig. S2E). MYC is a well-known oncogene22, and we found that in 63 pairs of patient tumor tissues (87.5%), MYC expression is more than 2 times higher than that in the corresponding native tissues. H3K27ac track on its enhancer was shown as an example (Fig. 1E). In the adjacent native tissues, MYCexpression was very low, and the H3K27ac signal on its enhancer and eRNA were close to the background; in tumor tissues, MYC was highly expressed, accompanied by elevation of H3K27ac on its enhancer and eRNA expression (Fig. 1E). One early study has reported active enhancers in multiple CRC cell lines and a few CRC samples33. Comparison of the two studies revealed that we identified 11796 different active enhancers in CRC, which was 32.4% of the total enhancers (Fig. 1F, G).

Identification of variant enhancer loci in tumor

To identify significant active enhancers specific in tumors, we first compared the enhancers of all samples and found that some tissues had a relatively low number of H3K27ac peaks (less than 2,500), or variant enhancer loci (VELs, less than 500) compared with the corresponding adjacent tissues (Supplementary Fig. S3A and B). We considered it was probably due to the sampling problem. Since when we were using random pieces of clinical dissected tissues for experiments, it was possible that some tissue samples only contained very few cells, or a large part of their cells actually were very similar as their corresponding ones. So, we ruled them out in the following statistical analysis. We totally identified 6690 significant VELs, including 5590 gain VELs and 1100 lost VELs (Fig. 2A, Supplementary Data 3–5) and the pipeline was shown in Supplementary Fig. S3C. At the recurrence threshold of 14 and 19 patients, 95% of the gain and lost VELs achieved statistical significance (q-value<0.1, paired t-test, with Benjamini-Hochberg correction; Fig. 2B). Supporting the reliability of these analysis, significant gain VELs exhibited higher H3K27ac level in tumors than native tissues, and opposite in significant lost VELs (Supplementary Fig. S3D and E). Meanwhile, genes associated with gain VELs showed elevated expression in tumors in comparison with native tissues, while lost VEL-associated genes were broadly repressed in tumors; and the magnitude of the change in expression positively correlated with the number of VELs per gene (Supplementary Fig. S3F). The gain VELs close to IL20RA and FOXQ1 and the lost VELs close to PPARGC1B were shown as representative (Fig. 2C, Supplementary Fig. S3G and H). The identified VELs could nicely distinguish the native and tumor tissues (Fig. 2A). Human disease ontology and GO analysis showed that the associated genes of gain VELs were highly related with CRC (Fig. 2D & Supplementary Fig. S3I), while those of lost VELs were related with normal colon functions (Supplementary Fig. S3J). To further evaluate the potential of H3K27ac or enhancer information in distinguishing tumor and normal tissues, we clustered the adjacent native and tumor tissues with PCA using the information of gene expression, H3K27ac on enhancers, and H3K4me3 on promoters. The adjacent and tumor samples were nicely distinguished using both gene expression or significant enhancers, while the different H3K4me3 peaks did not (Fig. 2E–G), suggesting enhancer information is useful for tumor identification.

A Relative H3K27ac signals of lost and gain VELs in all tumor and native tissues. B The required recurrence for gain and lost VELs meeting statistical significance (p.adj <0.05). The two vertical dashed lines at left highlights the recurrence of gain and lost VELs when achieve the cut-off (0.95, black dashed line) of significant percentage, and the two lines at right highlights the highest recurrence in tumor or native tissue of gain and lost VELs, respectively. For gain VELs, when recurrence reach to 14, the percentage of significant VELs is 95.635%; for lost VELs, when recurrence reaches to 19, the percentage of significant VELs is 96.273. p.adj indicates the BH adjusted t-test p-value. A two-sided test was utilized here. C Representative H3K27ac tracks of gain VEL on IL20RAloci. D The human disease ontology in which gain VELs participated detected by GREAT software (version 3.0.0). The red bars represent CRC-related diseases and the black bars represent other diseases. E–G PCA analyses to classify tumor and native tissues using gene expression (E), all significant enhancers (F), and promoters (G) information identified using our patient RNA-seq and ChIP-seq data.

Full size image

The enhancer features of CRC subgroups

To further investigate the enhancer features in CRC subgroups, we utilized one of the common approaches, the consensus molecular subtypes (CMS) classification of CRC tumors35, and classified patients into four subgroups (Fig. 3A and B, Supplementary Fig. S4B–C). The correlation analysis based on the identified VELs and subgroup-specific VELs showed that the tissues of the CMS2 group had the highest correlation, while CMS4 was the lowest (Fig. 3C and D, Supplementary Fig. S5), suggesting CMS4 might be more heterogeneous than others. Interestingly, when comparing the enhancers among the four subgroups, we found that CMS2 had the largest number of active enhancers, significant gain VELs, and specific gain VELs (Fig. 3E–G, Supplementary Fig. S6A). The average H3K27ac signal of gain VELs in CMS2 was also higher than the other three (Fig. 3F). While, the four groups had no big difference at the amount of gain VELs among individual samples, as well as the correlation of H3K27ac across the genome (Supplementary Fig. S6B–D). The above study indicates that the CMS2 group is more homogenous than others; and it has more specific active enhancers, which might be a novel feature for it.

A The patient identifier of members in four CMS subgroups. B The consensus molecular subtypes (CMS) classification of CRC samples using R package CMScaller. C Correlation of H3K27ac signal on the regions of gain VELs in all tumor samples of CMS1–4 subgroups. Correlations were calculated by the Spearman correlation coefficient. D The required recurrence for gain VELs in each CMS subgroup to meet statistical significance (p.adj <0.05) at different cut-offs. The dashed lines highlight the recurrence of gain VELs when achieve the cut-off (0.9, black dashed line) of a significant percentage. When recurrence reach to the cut-off, the significant percentage of CMS1 = 93.942%, CMS2 = 93.615%, CMS3 = 95.741% and CMS4 = 93.548%. p.adj indicates the BH adjusted t-test p-value. A two-sided test was utilized here. E The number of subgroup significant gain VELs in four CMS subgroups. F The average H3K27ac signal (RPM) on the regions of gain VELs in four CMS subgroups. G The number of subgroup-specific gain VELs in each CMS. The subgroup-specific gain VELs were identified when the mean RPM of one VEL in one CMS subgroup was 1.5 times higher than the other three. H Functional annotation of target genes associated with CMS2 specific gain VELs based on their significant overlap with gene sets annotated in Gene Ontology (Biological Process) and pathway database (Reactome). I Meta tracks of normalized H3K27ac on CEL and DPEP1 gene loci in four CMS subgroups.

Full size image

Then, we studied the function of VEL-associated genes, and identified enhancers and genes specifically activated in each subgroup. Some representative enhancers and genes for each group were shown (Supplementary Fig. S6E–L). For CMS2, we found its specific gain-VEL-associated genes are mainly involved in WNT signaling, cell migration, and lipid metabolic process (Fig. 3H, Supplementary Fig. S6M). Activation of WNT signaling and enhanced cell migration are expected, since APC is one of the most frequent mutated genes in CRC and cell migration is a hallmark for cancer cells36,37. Lipid metabolism was linked with CRC but ambiguous results from different groups exist38,39,40. Our analysis suggested dysregulation of lipid metabolic homeostasis is possibly associated with certain CRC subgroups. Some VELs of the CMS2 subgroup and their associated genes were shown, including those involved in lipid metabolism, such as CEL and DPEP1 (Fig. 3I, Supplementary Fig. S7).

Analysis and verification of variant super enhancer loci

Activation of oncogene-associated super-enhancers is one of the important features for cancer6. Using similar approaches as VEL identification, we identified the variant super-enhancers loci (VSEL) in tumor tissues, including 334 gain VSELs and 121 lost VSELs, among which several well-known oncogenic targets were identified, such as MYC, VEGFA, and LIF(Fig. 4A, Supplementary Fig. S8A and B, Supplementary Data 6&7). H3K27ac level on the gain VSELs were significantly increased and decreased on the lost VSELs as expected (Fig. 4B and C). As expected, genes associated with gain VSELs expressed higher than those associated with lost VSELs (Supplementary Fig. S8C). We utilized the H3K27ac value on VSELs to cluster CRC patients, together with normal intestinal tissues (Supplementary Fig. S8D and E). The analysis distinguished the native and tumor tissues, and classified CRC patients into three subgroups (Supplementary Fig. S8D and E). The results suggest VSELs might be useful for CRC classification. We compared the two classification approaches, and found that most of the CMS2 group samples were classified into G1 subgroup based on VSELs (Supplementary Fig. S8F), supporting our previous conclusion that the CMS2 group is more homogenous than others.

A The genes associated with top super-enhancers (SEs) ranked by recurrence. Red dots represent tumor-specific SE genes and blue dots represent native tissue-specific SE genes. Top 10 tumor and native tissue-specific genes were listed. B–C The average H3K27ac signal (RPM) at the regions of gain VSELs (B) and lost VSELs (C) in tumor and native tissues. D Meta normalized H3K27ac tracks at IER3 gene loci. The green track on the top represents the H3K27ac signal in HCT116, and the black and grey lines at the bottom represent the average signal of tumor and native tissues, respectively. The pink lines indicate the target positions of dCas9-KRAB sgRNAs. E Bar plot showing the relative mRNA level of LIF, SLC7A5, CYP2S1, PHF19, RNF43, CEBPB, TBC1D16, TNFRSF6B, VEGFA, and IER3 in control and sgRNA groups (n = 3). A sgRNA control targeting EGFP was used as a control in the following experiments. *p < 0.05. F Transwell assays for HCT116 cell lines stably transfected with dCas9-KRAB sgRNAs of the enhancers mentioned in Fig. 5E (n = 3). G–I Xenograft experiments in nude mice were performed with HCT116 stable cells expressing the indicated sgRNAs. The tumors were pictured (G), and their volume and weight were shown (H&I). n = 9 for all groups. Data are presented as mean values ± SEM. Statistical analysis was performed using a two-sided Student t test. p value was labelled on the corresponding items.

Full size image

To experimentally verify the functions of the identified VSELs, we compared the H3K27ac profiles on top gain VSELs of CRC tissues with those in HCT116 cells. The gain VSELs appearing in HCT116 were chosen and the dCas9-KRAB system was utilized to repress enhancer activity41. Since one SE usually covers a relatively large chromatin region, we design multiple sgRNAs for each SE and made stable cell lines in HCT116, RKO or SW620, three CRC cell lines (Fig. 4D, Supplementary Fig. S9, Supplementary Data 8). Totally we analyzed 11 SEs, among which 10 SEs were found to regulate the expression of their proximal genes in multiple cell lines, including IER3, LIF, SLC7A5, CYP2S1, PHF19, RNF43, CEBPB, TBC1D16, TNFRSF6B, and VEGFA (Fig. 4A and E, Supplementary Fig. S9). For some enhancers, the expression of multiple close genes was repressed (Sup. Fig. S9). We also measured the H3K27ac level on some loci and found that dCas9-KRAB/sgRNA effectively repressed H3K27ac on the enhancers of CEBPE, CYP2S1, IER3, PHF19, RNF43, and TBC1D16(Supplementary Fig. S10A). In comparison with their expression in CRC patient tissues Supplementary Fig. S10B), our data indicated that the above SEs regulate the expression of the corresponding genes. Moreover, we established stable cell lines of repressed enhancers in HCT116 and studied their proliferation and migration ability. The difference of proliferation and cell cycle was not very significant (Supplementary Fig. S11A and B), however, quite a few cell lines exhibited attenuated migration ability in multiple cell lines, including PHF19, LIF, SLC7A5, CYP2S1, RNF43, VEGFA, and TBC1D16 (Fig. 4F, Supplementary Fig. S11C–F). To further investigate the functions of the above SEs in CRC, we performed xenograft assays with stable HCT116 cells, which showed that repression of PHF19 and TBC1D16 SEs significantly reduced tumor growth, and the expression of target genes was confirmed in randomly selected tumors (Fig. 4G–I, Supplementary Fig. S11G–K). When SLC7A5 SE was repressed, tumor growth also showed a trend of reduction, but the difference was not significant enough (Fig. 4G–I).

Predication and verification of functional transcription factors

To investigate the potential transcription factors (TFs) playing key roles in CRC, the DNA sequences of VELs was used for prediction by HOMER software. The top hits of gain and lost VELs were listed (Fig. 5A and Supplementary Fig. S12A and B). The hypothesis of core regulatory circuitry was raised to identify core TFs in cells42,43. To improve TF prediction, we utilized the method to identify key TFs in CRC tissues (Fig. 5B&C, Supplementary Fig. S12C). ASCL2 was predicted as a CRC-specific TF with the highest score (Fig. 5C). The H3K27ac level of ASCL2 enhancer greatly increased in tumors (Supplementary Fig. S12D), and the gene expression analysis based on TCGA datasets suggested ASCL2 was highly expressed specific in colorectal cancer (Supplementary Fig. S12E). These suggest ASCL2 is a key TF in CRC, which is consistent with the previous publications44,45,46. Since ASCL2 has been characterized in CRC, we did not further investigate its roles.

A DNA motifs enriched within nucleosome-free regions (NFRs) of tumor gain VELs determined by HOMER motif analysis. B Heatmap of transcription factors ranked by predicted core regulatory circuitry (CRC) total degrees (Tumor - Native tissue). Top 30 tumor and native-specific TFs were listed. C Scatter plot showing the total degree (Tumor-Native tissue) and expression FC (Tumor/Native tissue) of the specific TFs listed in Fig. 5B. Blue dots represent top 30 tumor specific TFs, and red dots represent the top 30 native-specific TFs. Circle size indicates the mean expression (FPKM) of TFs in its specific tissue. D, E Transwell assay for HCT116 cell line with KLF3 knockdown. KLF3 was measured with western blotting. p-value = 1.06E−4 (sgKLF3.1), and 1.55E−4 (sgKLF3.2). n = 3. Data are presented as mean values ± SEM. F–H KLF3 stably knockdown HCT116 cells by sgRNA were injected into nude mice (10^6 cell pear mouse, n = 10). Tumors were pictured (F), and tumor growth curve (G) and weight (H) were shown as mean (± SEM). P-value in G for sgKLF3.1 and KLF3.2 are 0.0261 and 0.0296, in H for sgKLF3.1 and KLF3.2 are 0.0256 and 0.0355, respectively. Statistical analysis was performed using a two-sided Student t test. *p < 0.05, **p < 0.01.

Full size image

Combining the above results and the published literatures, we selected around 10 TFs for experimental verification. Unfortunately, some genes express low in the selected cell lines or knockdown experiments did not work for some genes. Eventually, we only got 4 genes successfully knocked down, including KLF3 and MAFK, two novel TFs, MAZ and RUNX1, two recently reported TFs functioning in CRC but not well characterized46,47,48,49. The significance of the selected TFs in each patient was calculated (Supplementary Fig. S12F). Knockdown of these genes did not affect cell proliferation in cell proliferation assay; while KLF3, MAZ, and RUNX1 knockdown by siRNAs, or KLF3 knockdown by CRISPR/sgRNAs, repressed cell migration, but not MAFK (Fig. 5D and E, Supplementary Fig. S13, S14A and B). Our results identified KLF3 as a TF involved in CRC, and confirmed the reported roles of RUNX1 and MAZ. MAFK was found at the side of native tissues in core regulatory circuitry analysis (Fig. 5B), and our experimental results did not find its role in proliferation or migration, suggesting MAFK is not involved in CRC although it was predicted in DNA motif assay. To further confirm the function of KLF3, we knocked it down with CRISPR/sgRNAs in HCT15 and RKO cell lines, and the similar results were observed (Supplementary Fig S14C–F). Xenograft assay was then performed with KLF3-deficient HCT116, and the result indicated that KLF3 knockdown repressed tumor growth in mice (Fig. 5F–H).

CRC is one of the most common cancers in the world. Though early screening greatly improves the curative ratio, novel classification approaches and drugs are still urged to be developed. The current study provides a comprehensive map of H3K27ac and active enhancers in CRC patients. The early studies used cell lines to determine CRC-specific enhancers33. Our study used paired patient tissues, which provided much more reliable data and identified many VELs and VSELs. Moreover, we experimentally confirmed the roles of more than 10 SEs in CRC. These provide important information for CRC research.

Our analysis predicted many TFs functioning in CRC, such as KLF3. ASCL2 was reported to be an oncogene in CRC, and our results showed that ASCL2 is highly expressed and H3K27ac on it increased in CRC tissues. Our results were consistent with published studies. We also experimentally confirmed the roles RUNX1 and MAZ, and identified KLF3 as an oncogenic TF in CRC. Further studies about these TFs will provide important information for CRC research. Loss of KLF3 was previously shown to be associated with more aggressive CRC50. We then analyzed the expression of KLF3 in our data and TCGA data. We did not observe a significant difference between the adjacent and cancer tissues in our collection; among TCGA samples, KLF3 expression was higher in tumor tissues but not significant enough. We then compared the expression between the adjacent and cancer tissues of four CMS groups. Interestingly, KLF3 is lower expressed in the cancer tissues of the CMS1 group, and higher expressed in the CMS3 group (Supplementary Fig. S14I). Although our sample number is a little bit small when split into four groups, the results suggest that KLF3 may play different roles in the four groups. It will be interesting to determine KLF3s function in each CRC subgroup.

Interestingly, verification of the identified VSELs and TFs indicated that most of them were more related with cell migration, but less with proliferation. It is possible that the normal intestinal epithelial cells are growing rapidly, and the change of epigenomic marks on migration-related genes is much bigger than proliferation genes, which caused the top candidates are mostly related to migration. Other possibilities could be that the CRC samples we collected are mostly at relative late stages (Supplementary Fig. S1A & Supplementary Data 1), and the tumor cells were probably at the metastasis stage or ready for it. The difference of enhancers and genes governing migration was probably much more significant than other genes between paired tissues, and the chosen VSELs were all among the most significant ones.

Interestingly, our study showed that SE repression sometimes caused down-regulation of multiple proximal genes, and the affected genes might vary in different cell lines (Supplementary Fig. S9F, G and I). It showed the complexity of transcription regulation by SEs, and indicated one fact that one SE is consisted of multiple TF binding sites and could be bound by multiple TFs, which may govern the transcription of multiple associated genes.

It was believed that the gain of H3K27ac on the oncogenic enhancers is a common feature for cancers21. Our analysis indicated that only the CMS2 group has the obvious feature of enhancer activity elevation, and the other subgroups have much less gained VELs. These indicate that in CRC, the global increase of active enhancers is an important feature for just one subgroup, not for all.

The homeostasis of lipid metabolism has been linked with CRC for many years, however, the detailed mechanisms is not clear and the use of statin analogues failed in CRC treatment38. We found that H3K27ac significantly increased on the enhancers of genes related with lipid metabolism. The bioinformatics analysis also pointed out that the gain VELs of the CMS2 subgroup were enriched with genes involved in lipid metabolic processes. So, it is possible that the dysregulation of lipid homeostasis is only associated with the CMS2 group, which should be explored by the future studies.