DNA methylation at enhancers identifies distinct breast cancer lineages

Abstract

Breast cancers exhibit genome-wide aberrant DNA methylation patterns. To investigate how these affect the transcriptome and which changes are linked to transformation or progression, we apply genome-wide expression-methylation quantitative trait loci (emQTL) analysis between DNA methylation and gene expression. On a whole genome scale, in cis and in trans, DNA methylation and gene expression have remarkably and reproducibly conserved patterns of association in three breast cancer cohorts (n = 104, n = 253 and n = 277). The expression-methylation quantitative trait loci associations form two main clusters; one relates to tumor infiltrating immune cell signatures and the other to estrogen receptor signaling. In the estrogen related cluster, using ChromHMM segmentation and transcription factor chromatin immunoprecipitation sequencing data, we identify transcriptional networks regulated in a cell lineage-specific manner by DNA methylation at enhancers. These networks are strongly dominated by ERα, FOXA1 or GATA3 and their targets were functionally validated using knockdown by small interfering RNA or GRO-seq analysis after transcriptional stimulation with estrogen.

Fig. 1

Identification of expression-methylation QTL (emQTL) a Unsupervised clustering of (−log) p-values of the emQTL by Pearson’s correlation and average linkage revealed two main clusters of CpG-gene pairs. Rows represent CpGs and columns represents genes. Yellow and light yellow spots show highly significant associations between CpG methylation and gene expression. b Density plot showing the degree of absolute co-expression of genes and co-methylation of CpGs in Cluster 1 (red) and Cluster 2 (blue). c, d Gene set enrichment analysis in Cluster 1 (c, n = 160) and Cluster 2 (d, n = 270) using MSigDB (H and C5 databases). The height of the bars represents the level of enrichment measured as a ratio between the number of genes overlapping an MSigDB H or C5 gene set over the expected frequency if such overlaps were to occur at random

Fig. 2

Genomic location of emQTL-CpGs according to ChromHMM and TF binding regions a Bar plot showing the enrichment of emQTL, Cluster 1 and Cluster 2 CpGs over the expected frequency of CpGs from the Illumina HumanMethylation450 across functional/regulatory regions of the genome as determined by MCF7 ChromHMM annotation. The height of the bars represents the level of enrichment measured as a ratio between the frequency of all emQTL-CpGs (white), CpGs in Cluster 1 (red) and CpGs in Cluster 2 (blue) overlapping a functional element over the expected frequency if such overlaps were to occur at random. Statistically significant enrichments (p < 0.05; hypergeometric test) are marked with an asterisk. b Enrichment analysis of CpGs in Cluster 2 across ERα (white), GATA3 (blue), FOXA1 (red) and CTCF (green) binding regions as determined by ChIP-seq analysis. Enrichment is calculated for different genomic regions determined by MCF7 ChromHMM annotation. For this analysis some ChromHMM annotations were collapsed into one as follow: Enhancer = ‘Enhancer’ and ‘Enhancer + CTCF’ and Promoter = ‘Promoter’, ‘Promoter + CTCF’ and ‘Poised Promoter’. The height of the bars represents the level of enrichment measured as a ratio between the frequencies of CpGs in Cluster 2 in each ChIP-seq peak at specific regulatory region over the expected frequency if such overlaps were to occur at random. Statistically significant (p < 0.05; hypergeometric test) enrichments are marked with an asterisk

Fig. 3

Differential DNA methylation of Cluster 2 CpGs in TF binding regions. Unsupervised clustering of DNA methylation levels of the CpGs in Cluster 2 for (a) the TCGA cohort (N = 609) and (b) the OSL2 cohort (N = 272). Annotation of the rows of the heatmap shows whether a CpG is located in ERα (yellow), FOXA1 (pink) or GATA3 (light blue) binding regions according to ChIP-seq experiments in MCF7 cells. Annotations of the column of the heatmap indicate histopathological features of the patients: PAM50 subtype, ER and PR status. c–d Average DNA methylation of c the 2097 CpGs in Cluster 2A and D) 1504 CpGs in Cluster 2B for samples in the TCGA cohort. Boxplots represent the average DNA methylation of these CpGs in ER positive (blue; n = 418), ER-negative tumors (red; n = 124) and adjacent normal tissue (green; n = 97). Kruskal-Wallis test p-values are denoted

Fig. 4

Differential expression of genes in Cluster 2. Unsupervised clustering of expression of the genes in Cluster 2 for (a) the OSL2 cohort (N = 272) and (b) the TCGA cohort (N = 528). Genes in rows are annotated if they are locally paired with CpG in Cluster 2 A (pink) or paired with CpG of Cluster 2B (light blue). Two main sub-clusters of genes are identified in the heatmap in a, and the genes in the heatmap in b are annotated correspondingly. Annotations of the column of the heatmap indicate histopathological features of the patients: PAM50 subtype, ER and PR status. c–d Average expression of genes in Cluster 2 locally paired with CpGs in Cluster 2 A (n = 58) and Cluster 2B (n = 31). A gene was considered to be locally paired with a CpG if it was situated not more than 10 kb from this CpG. Average expression in ER positive tumors (blue; n = 406), ER-negative tumors (red; n = 117) and adjacent normal tissue (green; n = 61) from the TCGA cohort. Kruskal–Wallis test p-values are denoted

Fig. 5

DNA methylation of ERα, FOXA1 and GATA3 binding regions and expression of their target genes in Cluster 2. a–c Average DNA methylation of CpGs in Cluster 2 and ERα a, FOXA1 b and GATA3 c binding regions defined by ChiP-seq peaks. Boxplots represent the average DNA methylation of these CpGs in ER positive (blue, n = 418), ER-negative tumors (red, n = 124) and adjacent normal tissue (green, n = 97). The average methylation of the Cluster2-CpGs in a TF binding site was significantly lower in ER positive patients compared to ER-negative and adjacent normal tissue. d–f Average gene expression of TF target genes in Cluster 2. d Estrogen (GRO-seq), e FOXA1 (siRNA) and f GATA3 (siRNA). Boxplots represent the average expression of the TF target genes in ER positive tumors (blue, n = 406), ER-negative tumors (red, n = 117) and adjacent normal tissue (green, n = 61). The average expression was significantly higher in ER positive tumors compared to ER-negative and adjacent normal tissue. Kruskal–Wallis test p-values are denoted

Fig. 6

Enhancer-promoter interaction and impact of TF binding on target gene expression. a Bar plot showing the enrichment of emQTL in ChIA-PET Pol2 loops for Cluster 1, Cluster 2, Cluster 2A and Cluster 2B. The height of the bars represents the level of enrichment measured as a ratio between the frequencies of emQTL (CpG–Gene pairs) found in the head and tail of Pol2 loops, over the expected frequency if such overlaps were to occur at random. Statistically significant enrichments (hypergeometric test, p-value < 0.05) are marked with an asterisk. b Example of overlap of emQTL (red arcs) and ChIA-PET Pol2 loops (blue arcs). Also shown are the location of ERα, FOXA1 and GATA3 binding regions. c, d dCas9 and ERα ChIP were performed in control (−gRNA and dCas9) or transfected MCF7 cells (gRNA E5 and dCas9), to assess the binding of each protein at enhancer or promoter. Statistically significant differences (t-test; two tails, p-value < 0.05) are marked with an asterisk. The data are presented as mean of three of independent replicates ± s.d. e mRNA levels of PGR were measured in control (-gRNA and dCas9) or transfected MCF7 cells (gRNA E5 and dCas9) by real-time PCR. Statistically significant differences (t-test; two tails, p-value < 0.05) are marked with an asterisk. The data are the mean of three of independent replicates ± s.d

Fig. 7

Circos plots showing the genomic location of all CpGs in emQTL with ESR1, FOXA1 and GATA3. Circos plot representing all associations between CpGs and ESR1 (a, d), FOXA1 (b, e), and GATA3 (c, f). Red lines represent negative associations a–c and blue lines represent positive correlations d–f. The outer ring indicates whether a CpGs is located in an enhancer determined by MCF7 ChromHMM annotation and the inner ring whether it is located in a binding region of the respective TF determined by ChIP-seq peaks

Fig. 8

Cluster 1 highlights a link between DNA methylation and lymphocyte infiltration. Unsupervised clustering of expression levels of the 160 genes in Cluster 1 from a the TCGA cohort (N = 528) and b the OSL2 cohort (N = 272). Annotations of the column indicate the level of lymphocytes infiltration, ER and PR status. Levels of lymphocyte infiltration were calculated from a set of genes expressed by lymphocyte characterized by the algorithm Nanodissect. c 68 tumor tissue samples were scored as low (n = 30), moderate (n = 22), high (n = 13) or severe (n = 3) inflammation by a pathologist based on the level of immune cell infiltration found in each tumor. Boxplot representing average expression of genes in Cluster 1 (n = 160) according to inflammation score. Kruskal–Wallis test p-value is denoted