Search

Additional file 1: Fig. S1. Flowchart describing data and the different steps of the analysis leading to the identification of 89,118 miRNA-methylation Quantitative Trait Loci (mimQTLs). Examples of negative and positive correlation between methylation at a CpG and expression of a miRNA are shown as scatterplots at the bottom. Fig. S2. Overview of the 89,118 miRNA-CpG associations found significant in both cohorts. The scatterplots show a) the –log10(Bonferroni-adjusted Spearman correlation p-values) of Oslo2 (x-axis) vs. TCGA (y-axis); b) Spearman correlation coefficients in Oslo2 (x-axis) vs. TCGA (y-axis). The histograms show the distribution of the correlation coefficients (Spearman’s rho) of all significant miRNA-CpG correlations in the Oslo2 (c) and TCGA (d) cohorts. Fig. S3. Number of associations and genomic positions of mimQTL miRNAs and CpGs. a) Barplot showing the number of negative (neg) and positive (pos) CpG correlations (cor) for the three different miRNA clusters. b) mimQTL Manhattan plot with genomic coordinates of CpGs (black or gray) and miRNAs (green) displayed along the x-axis, with the negative logarithm of the Bonferroni-corrected Spearman correlation p-value from Oslo2 on the y-axis. Each dot on the plot signifies a CpG or miRNA (CpGs are shown in two colors to distinguish the chromosomes more clearly). c) In cis mimQTL Manhattan plot displaying the chromosomal location (using the position of the CpG) along the x-axis of the 5125 mimQTLs found on the same chromosome (in cis). Each dot represents one mimQTL which is color-coded according to negative (black) or positive (green) miRNA-CpG correlation. The y-axis displays the negative logarithm of the Bonferroni-corrected Spearman correlation p-value from Oslo2. Fig. S4. Barplots showing the number of associations per miRNA or CpG. a) Barplot showing the number of CpG associations per miRNA (n = 119). Note that the y-axis is on log scale. b) Barplot showing the number of miRNA associations per CpG (n = 26,746). Fig. S5. Density plots showing the degree of CpG co-methylation or miRNA co-expression between cluster members (see Fig. 1 and Additional file 3 a, b) calculated by Spearman correlation. a) Correlation between CpG cluster members in the Oslo2 data. b) Correlation between CpG cluster members in the TCGA data. c) Correlation between miRNA cluster members in the Oslo2 data. d) Correlation between miRNA cluster members in the TCGA data. The dotted lines represent density plots of corresponding correlations expected by chance, i.e. correlations observed after randomly permuting the same data before performing correlation analyses. Fig. S6. Density plots showing the distribution of Spearman correlation coefficients between miRNA expression and selected variables for members of each of the miRNA clusters. a, b) miRNA expression-immune infiltration score [34] correlations for the Oslo2 (a) and TCGA (b) cohorts. c, d) miRNA expression-fibroblast infiltration score [35] correlations for the Oslo2 (c) and TCGA (d) cohorts. e, f) miRNA expression-ESR1 mRNA expression correlations for the Oslo2 (e) and TCGA (f) cohorts. Fig. S7. Heatmaps showing hierarchical clustering of miRNA expression levels (rows) from tumors (columns) of the Oslo2 (top) and TCGA (bottom) cohort. Clustering was performed using Euclidean distance and average linkage. Tumors are annotated with the following clinical/molecular classifications: PAM50 molecular subtypes (Luminal A (LumA), Luminal B (LumB), Basal-like (Basal), HER2-enriched (Her2), Normal-like (Normal); Lymphocyte infiltration (LI) group where tumors were divided into quartiles: 1 (low) – 4 (high); Fibroblast infiltration group (Fibro) where tumors were divided into quartiles: 1 (low) – 4 (high); Human epidermal growth factor receptor 2 (HER2) status; Estrogen receptor (ER) status. a, b) Clustering of miRNA cluster A expression (n = 23); c, d) Clustering of miRNA cluster B expression (n = 59); e, f) Clustering of miRNA cluster C expression (n = 37). Fig. S8. Heatmaps showing hierarchical clustering of methylation levels of CpG cluster 1 (a; n = 14,040) and CpG cluster 2 (b; n = 12,706) in the TCGA cohort (CpGs in rows and tumors in columns). Clustering was performed using Euclidean distance and average linkage. Tumors are annotated with the following clinical/molecular classifications: PAM50 molecular subtypes (Luminal A (LumA), Luminal B (LumB), Basal-like (Basal), HER2-enriched (Her2), Normal-like (Normal); Lymphocyte infiltration (LI) group where tumors were divided into quartiles: 1 (low) – 4 (high); Human epidermal growth factor receptor 2 (HER2) status; Estrogen receptor (ER) status. The CpGs are annotated according to overlap with regions annotated as “active intergenic enhancer” from ChromHMM of subtype-specific cell lines [37] with corresponding subtype colors. Fig. S9. Boxplot showing average DNA methylation of CpGs from cluster 1 in PAM50 subtypes of the Oslo2 cohort (Luminal A (LumA), Luminal B (LumB), Basal-like (Basal), HER2-enriched (Her2)). b) Boxplot showing average DNA methylation of CpGs from cluster 2 in the TCGA cohort when tumors were separated into quartile lymphocyte infiltration groups from low (1) to high (4) infiltration. c) Boxplot showing average DNA methylation of CpGs from cluster 2 in normal breast tissue (reduction mammoplasty, n = 17) or estrogen receptor (ER) positive (pos) or negative (neg) tumors of the Oslo2 cohort. d) Boxplot showing average DNA methylation of CpGs from cluster 2 in normal breast tissue (normal adjacent breast tissue, n = 97) or ER positive or negative tumors of the TCGA cohort. P-values resulting from Kruskal-Wallis tests indicated. Fig. S10. Boxplot showing DNA methylation of the hub CpG of miRNA cluster A (cg14270581; y-axis)) in ER positive (pos) and negative (neg) breast cancer cell lines and from different immune cell types (x-axis); B-cells, leukocytes (leuko), monocytes (mono) and T-cells. P-value resulting from Wilcoxon rank-sum test between cancer cell lines vs. immune cells is indicated. Fig. S11. Density plot showing the distribution of the Global Methylation Alteration (GMA) score in normal adjacent breast tissue (green), tumors (black) and tumors separated into estrogen receptor (ER) positive (pos) and negative (neg). Data from TCGA. Fig. S12. Top panel: Scatterplots showing on the x-axis mRNA expression of DNMT3A (left), DNMT3B (middle) and DNMT1 (right) vs. hsa-miR-29c-5p expression (y-axis) measured in 377 samples of the Oslo2 cohort. Bottom panel: Scatterplots showing on the x-axis protein expression of DNMT3A (left), DNMT3B (middle) and DNMT1 (right) vs. hsa-miR-29c-5p expression (y-axis) measured in 45 samples of the Oslo2 cohort. Each dot represents a tumor color-coded according to PAM50 subtype (Luminal A (LumA): dark blue; Luminal B (LumB): light blue; Basal-like (Basal): red; HER2-enriched (Her2): pink; Normal-like: green). Spearman correlation coefficient and p-value are indicated for each plot.

Breast cancer is the most prevalent type of cancer affecting females in world today and poses a serious public health problem worldwide. Even though the overall survival has increased significantly the last decades, the high incidence of breast cancer signifies the importance of improvements in diagnostics and treatment. Like all types of cancers, breast cancer is a result of accumulation of genetic and epigenetic alterations that leads to repression of tumor suppressor genes and activation of oncogenes. Over the past decades, several studies have highlighted alterations in DNA methylation patterns as hallmark events in many cancer types including breast cancer. Among the major types of breast cancers, the estrogen receptor positive tumors which accounts for 70 % of all breast cancers tend to display more pronounced changes in their DNA methylation landscape than the ER negative tumors when compared to normal adjacent tissue. It is well known that alterations in DNA methylation may affect the expression of genes, depending on where the changes occur. For instance, changes in DNA methylation at CpGs in cis-regulatory regions such as promotors tend to repress the expression of its associated gene. Furthermore, CpGs as far as 100 kb away from the transcription start site have been demonstrated to be associated with gene expression. Therefore, CpGs in intergenic and enhancer regions may play key roles in breast cancer pathogenesis through the regulation of expression of their associated genes. Enhancer methylation has been shown experimentally to be associated with gene expression, and these genomic regions are considered to be the most differentially methylated genomic regions during carcinogenesis and cancer progression. Enhancers are known to bind cell type specific proteins called transcription factors (TFs) which are proteins involved in the regulation of gene transcription. However, the role of DNA methylation at enhancer regions regarding TF binding and breast cancer pathogenesis is still not fully understood. Genome-wide expression-methylation quantitative trait loci (emQTL) analysis have previously been shown to identify significant correlations between the level of DNA methylation at CpG sites and gene expression due to intertumoral heterogeneity within ER positive and ER negative tumors. It has also been shown to be a valuable tool in the identification of key gene regulatory networks involved in breast cancer pathogenesis. To take this further, the same approach was applied to the ER positive breast tumors only, to investigate whether any differences within the ER positive tumors in respect to DNA methylation and gene expression could be observed. The study resulted in the identification of CpG-gene pairs in which the level of DNA methylation was significantly correlated with gene expression. Hierarchical clustering of the significant associations led to the discovery of a previously undiscovered cluster of CpG-gene associations. Gene set enrichment analysis indicated an enrichment of the genes in EMT-related processes, while the CpGs were highly enriched in enhancer regions. The CpGs in this EMT-cluster was divided into CpG-cluster A and CpG-cluster B based on whether their mean methylation value was more or less than 0.5 respectively. The CpGs in both clusters were shown to be differentially methylated among the ER positive tumors. Further characterization of the CpG-clusters by ChIP-seq peaks enrichment analysis revealed that CpG-cluster A CpGs were enriched within ChIP-seq peaks of TFs associated with EMT such as TEAD1, FOSL1, TWIST1, SIX2, YAP1 and PPARG. To investigate whether the difference in DNA methylation was associated with any phenotypic feature associated with EMT in the tumors, an EMT score was utilized and correlated with the mean methylation of CpG-cluster A CpGs. The mean methylation was negatively correlated with the EMT score, meaning that lower methylation was associated with more mesenchymal-like characteristic of the tumor samples. These findings suggest that EMT-related CpG-gene pairs discovered in this study are associated with gene regulatory networks wired by EMT related TFs through a relationship between DNA methylation at their target DNA binding regions in enhancers, and gene expression of their target genes. This study highlights the CpGs identified as potential contributors to EMT-related cancer pathogenesis in the ER positive breast tumors and constitute interesting regions for further investigations. However, these in silico findings still require better validation in vitro. The identification of cancer-causing epigenetic changes may open up possibilities of targeted treatment by utilization of technologies such as CRISPR to edit epigenetic cancer-causing mutations to inhibit tumor growth in the future.