Epigenetic alterations at distal enhancers are linked to proliferation in human breast cancer

Abstract

Aberrant DNA methylation is an early event in breast carcinogenesis and plays a critical role in regulating gene expression. Here, we perform genome-wide expression-methylation Quantitative Trait Loci (emQTL) analysis through the integration of DNA methylation and gene expression to identify disease-driving pathways under epigenetic control. By grouping the emQTLs using biclustering we identify associations representing important biological processes associated with breast cancer pathogenesis including regulation of proliferation and tumor-infiltrating fibroblasts. We report genome-wide loss of enhancer methylation at binding sites of proliferation-driving transcription factors including CEBP-β, FOSL1, and FOSL2 with concomitant high expression of proliferation-related genes in aggressive breast tumors as we confirm with scRNA-seq. The identified emQTL-CpGs and genes were found connected through chromatin loops, indicating that proliferation in breast tumors is under epigenetic regulation by DNA methylation. Interestingly, the associations between enhancer methylation and proliferation-related gene expression were also observed within known subtypes of breast cancer, suggesting a common role of epigenetic regulation of proliferation. Taken together, we show that proliferation in breast cancer is linked to loss of methylation at specific enhancers and transcription factor binding and gene activation through chromatin looping.

Figure 1.

Study overview. Flowchart showing the pipeline used for identification of CpG-gene associations (emQTLs) and methods used for emQTL grouping and characterization.

Figure 2.

Identification and characterization of the emQTL biclusters. (A) Line chart showing the average MSR score for the biclusters obtained by biclustering of the inverse correlation coefficients obtained from the emQTL analysis in OSL2 when the number of biclusters k were set to be a number between 2 and 20. (B) Heatmap showing the inverse correlation coefficients of the emQTL-CpGs (n = 44 263) and emQTL-genes (n = 4904) from OSL2 after biclustering. Rows represent CpGs and columns represent genes. Each of the five biclusters is annotated. Blue points indicate negative correlations between the variables while red points represent positive correlations. White points indicate little or no correlation. (C) GSEA of the genes in Bicluster 1 (n = 1085), Bicluster 2 (n = 870), Bicluster 3 (n = 936), Bicluster 4 (n = 1087) and Bicluster 5 (n = 926) using gene sets obtained from the MSigDB. The length of the bars represents the log-transformed Benjamini-Hochberg (BH) corrected P-values obtained by hypergeometric distribution. Red bars indicate Hallmark gene sets while GO biological process, GO molecular function, and GO cellular compartment GO gene sub-collections are colored in orange, green and blue, respectively. Overlap between the gene list of the bicluster and each MSigDB gene set is annotated within each bar.

Figure 3.

Functional characterization of the emQTL-CpGs in the cell cycle bicluster. (A) Bar plot the showing enrichment of the cell cycle bicluster-CpGs in ChromHMM-defined genomic regions by subtype. The length of the bars represents the log-transformed BH corrected P-values. The color gradient of the bars represents fold enrichment (FE) in which a red color indicates FE close to 3.5 while white bars are genomic regions by subtype with FE close to 0. An enrichment was considered to be significant if the BH-corrected P-value was less than 0.05. (B) UpSet plot showing the overlap between CpGs in the cell cycle bicluster found within ChromHMM-defined active intergenic enhancer regions by breast cancer subtype. (C) Bar plot representing enrichment of the cell cycle bicluster-CpGs in the binding region of specific TFs according to UniBind. Bar length displays the log-transformed BH-corrected P-value obtained by hypergeometric testing for each TF. Red color indicates FE close to 3.5 while a white color indicates FE close to 0. (D) Unsupervised hierarchical clustering of DNA methylation levels of the 8641 cell cycle bicluster-CpGs for the tumor tissue from OSL2 with PAM50 status available (n = 272). Rows represent CpGs and columns represent histopathological features including PAM50 subtype and ER status of the tumor samples. Red points indicate methylated CpGs while blue points represent unmethylated CpGs. Boxplots showing the average DNA methylation of the cell cycle bicluster-CpGs in (E) OSL2 (n = 272) and (F) TCGA (n = 562) by PAM50 subtype. Kruskal-Wallis test P-values are denoted in the lower-left corner.

Figure 4.

Expression of genes in the cell cycle bicluster. (A) Unsupervised clustering of the expression levels of the 1085 genes in the cell cycle bicluster for the tumors in the OSL2 breast cancer cohort (n = 272). Rows represent genes and columns represent samples annotated with histopathological features including PAM50 subtype and ER status. Red color indicates high expression levels and blue color indicates low. Boxplots showing the average expression of genes in the cell cycle bicluster in the (B) OSL2 (n = 272) and (C) TCGA (n = 981) breast cancer cohorts. Kruskal-Wallis test P-values are denoted. Scatterplots showing the association between average DNA methylation of the cell cycle bicluster-CpGs versus average expression of the genes contained within the same bicluster by ER status in the OSL2 (D) and TCGA (E) breast cancer cohorts. Pearson correlation coefficients and P-values are denoted and colored by ER status.

Figure 5.

DNA methylation at enhancers facilitates target gene expression through enhancer-promoter interactions. (A) Bar plot showing the enrichment of emQTLs in ChIA-PET Pol2 loops and IM-PET loops for the ER+ MCF7 and ER– HCC1954 breast cancer cell lines, respectively. Bar height represents the enrichment level measured as the ratio between the frequency of emQTLs (CpG-gene pairs) found in the head and tail of a loop over the expected frequency if such overlaps were to occur at random. Enrichments that are statistically significant (hypergeometric test, BH corrected P-value < 0.05) are marked with an asterisk. (B) Enhancer hypomethylation at specific enhancers allows TF binding and the transcriptional activation of enhancer target genes through physical enhancer-promoter interactions by chromatin looping. (C) An example of a potential proliferation-promoting alteration in which the CpG (cg00733115) has been found in one foot of a ChIA-PET Pol2 loop (red arc) and a gene associated with proliferation (PIM1) is found in the other. Annotations for active intergenic enhancer regions and active promoters according to ChromHMM that are conserved across the cell lines of a similar subtype are shown in green and blue color respectively by breast cancer subtype. The binding sites of FOS, FOSL1/2 are also shown. (D) Scatterplot showing the association between DNA methylation at the emQTL-CpG cg00733115 and its associated gene (PIM1) by ER status in OSL2. Pearson's correlation coefficients and P-values are denoted.

Figure 6.

Expression of genes in the cell cycle bicluster associates with prognosis. Kaplan-Meier survival curves for the cell cycle bicluster in METABRIC cohort, for Luminal A (A), Luminal B (B), Basal-like (C), Normal-like (D), Her2-enriched (E) and all breast cancer subtypes (F). Tumors were divided into two groups based on the median of the average expression of genes in the cell cycle bicluster. P-values obtained by log-rank test are denoted.

Figure 7.

The EMT bicluster highlights an association between DNA methylation and fibroblast infiltration. (A) Heatmap showing the unsupervised clustering of the expression levels of the 936 genes contained within the EMT bicluster for 272 tumor samples from the OSL2 cohort. Rows represent genes and columns represent tumor samples annotated by histopathological features including PAM50 subtype and ER status. The tumor samples were divided into quartile groups based on the severity of fibroblast infiltration according to the relative amount of fibroblast in the tumor samples estimated by xCell. Differences in expression of the EMT bicluster genes between the quartile groups are shown for the OSL2 (B) and TCGA (C) cohorts. Each quartile group consisted of 68 tumor samples in OSL2 and 139 samples in TCGA. Boxplots showing the average DNA methylation at the 6910 CpGs contained within the EMT bicluster according to fibroblast infiltration score in (D) OSL2 (n = 272) and (E) TCGA (n = 556). Average DNA methylation values for these CpGs for in the PMC42-LA before and after EGF induced EMT. Fibroblasts, and the ER+ MCF7 and ER– MDAMB436 breast cancer cell lines are also included. Kruskal–Wallis test P-values are denoted in the bottom left corner.

Figure 8.

Cell-type-specific expression of the emQTL-bicluster genes. Combined Uniform Manifold Approximation and Projection (UMAP) plot for all 14 breast cancer samples annotated by cell type (A). Dot plots showing the expression of the selected genes for each bicluster for B all patients and patients with a (C) Luminal A, (D) Luminal B, (E) Her2-enriched and (F) TNBC. The size of the dot depicts the percentage of cells within each class and the intensity of the color shows the average expression level of each class.