Abstract

1. Introduction

Breast cancer is the most frequently diagnosed cancer in women in the United States accounting for 26% of newly diagnosed cases in 2008. It is the major cause of death among adult women, and lifetime risk of dying from breast cancer is 33 per thousand among women from high income countries. Women in high income countries have a higher risk than women from middle or low income countries, reflecting different exposure to known risk factors such as hormonal exposure and weight, age at menarche, number of pregnancies brought to term, and extent of breast feeding. Only about 10% of breast cancer cases are attributed to known hereditary factors, such as mutations in BRCA1 and BRCA2. Breast carcinomas are classified using clinical (tumor size, lymph node status) and histo‐pathological (grade, hormone receptor status) metrics to evaluate likely aggressiveness and identify optimal courses of treatment. The classical prognostic factors that are typically used in the clinic are node status, tumor size and tumor grade. Estrogen receptor status and ErbB2 status are therapy response predictors but do not have significant prognostic value independent of therapy.

For some time, the major characteristic differentiating breast carcinomas for treatment purposes has been estrogen receptor (ER) status. More recently, sophisticated molecular analyses including measurements of gene expression and genomic aberrations have refined breast cancer classfication. Most classifications mirror the separation seen by ER status but with additional information identifying smaller groups with homogenous molecular alterations and/or clinical behavior. The most widely accepted molecular classification of breast carcinomas is the “Intrinsic Classification” based on gene expression first proposed by Perou and Sorlie in 2000 (Perou et al., 2000; Sorlie et al., 2003). They identified five major subtypes of breast tumors: Luminal A and Luminal B (dominated by ER positive samples and with expression of genes typical of glandular epithelium); Basal‐like (ER negative samples with expression of myoepithelial associated genes); ErbB2+ (heterogeneous ER status, but often amplified for ErbB2+ and nearby genes, and a strong resemblance to Basal‐like expression); and Normal‐like (tumors with expression patterns close to normal breast tissue). Of the five subclasses, the Luminal A and Basal‐like subtype are the most clearly defined. Several studies have shown that they represent tumors with different aberrations at the genomic level as well, and these groups correspond reasonably well to clinical characterization on the basis of ER and HER2 status, as well as proliferation markers or histological grade (Bergamaschi et al., 2006). The intrinsic subtypes have also been associated with different prognostic implications with the Luminal A subtype having better prognosis than the Basal‐like and ErbB2+ subtypes.

Several recent studies have reported on the epigenetic influences in breast cancer (Hinshelwood and Clark, 2008). These include several classically studied breast cancer genes such as ESR1, CDH1 and CDKN2A (Birgisdottir et al., 2006; Caldeira et al., 2006). Flanagan et al., (Flanagan, Cocciardi et al.) used Affymetrix promoter arrays and methyl DNA immune‐precipitation to study 33 familial breast cancers. They found DNA methylation patterns are significantly associated with BRCA mutation status, although they could not determine association with subtypes in their small sample set.

We set out to determine if the gene expression‐based subtypes have underlying epigenetic differences. Epigenetic modifications both at the chromatin and DNA level affect the structure and the expression of genes encoded in the DNA. The most widely studied epigenetic modification is the cytosine methylation in the context of the dinucleotide CpG. In embryonic stem cells such modifications is of major importance in regulating genes important for cell differentiation and function. Altered regulation of CpG methylation is also implicated in many diseases. Specifically, in cancer, methylation of CpG islands proximal to tumor suppressors such as p16, RASSF1A, and BRCA1, is a frequent event (Merlo et al., 1995; Rice et al., 1998; Dammann et al., 2000). Several high throughput microarray based methods have been described that employ methylation‐sensitive restriction enzymes for detection of methylation (Huang et al., 1999; Khulan et al., 2006; Ordway et al., 2006; Irizarry et al., 2008). We recently reported one such method, called Methylation Oligonucleotide Microarray Analysis (MOMA) (Kamalakaran et al., 2009). MOMA allows for high throughput analysis of thousands of genomic loci including most CpG islands, and requires as little as 100 ng of sample DNA. The large number loci that can be interrogated using MOMA allows an unbiased investigation of genome‐wide methylation patterns as well as an in‐depth search for loci whose methylation have clinical importance in breast cancer.

The Oslo Micrometastases Study (OMS) (Wiedswang et al., 2003) has been the subject of one of the largest concentration of parallel molecular analysis techniques on a clinical dataset. The full study comprises over 900 breast cancer cases with information about presence of disseminated tumor cells (DTC), and associated clinical information such as node status, tumor size, estrogen receptor status, grade as well as a median of 85 months in follow‐up. A subset of approximately 140 patients is represented with fresh frozen samples from the primary tumor, matched blood, and micrometastases and has been used in parallel pilot studies of immunohistochemisty (Bergamaschi et al., 2006), whole genome mRNA expression (Naume et al., 2007), arrayCGH, whole genome SNP and SNP‐CGH, whole genome miRNA expression analyses and high throughput sequencing. In this paper we report the results of one of two independent DNA methylation studies performed on this information‐rich set employing the genome‐wide MOMA method. The MOMA method provides a global view of methylation state of the genome, covers most of the CpG islands and provides resolution at a sub‐CpG island level dependant on restriction fragment generated by representation. However, some classically studied CpG islands associated with cancer‐related genes such as p16 and RASSF1A are missed by MOMA because of excessive fragmentation by the representational enzyme. To address this issue, a complementary and parallel study was carried out using Illumina Golden Gate array methodology and the results are reported in the concurrently submitted paper by Rønneberg et al The Illumina method addresses a much smaller but focused portion of the genome (807 cancer‐related genes), but provides single nucleotide resolution for 1505 CpG sites in these genes. Together, these two studies provide a detailed survey of epigenetics in breast cancer, their relationship to the molecular subtypes and other clinical factors such as hormone status and TP53 mutational status and their implications for relapse risk.

2. Results

We performed MOMA (Methylation Oligonucleotide Microarray Analysis) on 119 breast samples (108 frozen breast tumors plus 11 normal adjacent samples collected during surgery). The primary tumors from the Oslo Micrometastases Study (OMS) (Wiedswang et al., 2003) were from stage I to stage III patients, and came with a variety of clinical, pathological and molecular data and have been described previously. A summary of the clinical information on the samples is provided in Table 1. MOMA has been described and validated using cell lines, breast and ovarian tumor tissues previously (Kamalakaran et al., 2009). Briefly, the sample DNA is digested with MspI to make a representation, ligated with adapters and split into two pools. One pool is restricted with the methylation‐specific restriction enzyme McrBC, while the other pool is a mock treated control. The two pools are then amplified, fluorescently labelled and hybridized onto a microarray. The ratio of intensity of the mock‐treated sample over McrBC‐treated sample provides an estimate of methylation. The methylation states of each fragment are determined by using an expectation maximization algorithm to assign probability of each fragment belonging into each of three distinct states – unmethylated (−1 state), partially methylated (0 state), and methylated (+1 state). A fragment is determined to be differentially methylated if it switches states from one state to another across samples (−1 to 0 or 0 to +1 or −1 to +1).

Table 1

Clinical Characteristics of tumors profiled in this study.

Sample Category	Sample Number
Samples
1. Tumors	108
2. Normal	11
Expression Subtype
1. Luminal A	40
2. Luminal B	13
3. ErbB2‐Like	19
4. Basal	12
5. Normal‐Like	14
6. Not Available	10
TOTAL	108
Hormone Receptor Status
1. Positive	66
2. Negative	35
3. Not Available	7
TOTAL	108
Tumor Grade
1. Grade I	14
2. Grade II	52
3. Grade III	40
4. Not Available	2
TOTAL	108
Tumor Size (Category)
1. pT1	46
2. pT2	49
3. pT3	6
4. pT4	4
5. Not Available	3
TOTAL	108
Lymph Node Status
1. pN0	43
2. pN1	31
3. pN2	18
4. pN3	8
5. Not Available	8
6. TOTAL	108

Open in a separate window

We first identified loci whose methylation state differed in breast tumors when compared to normal tissue to determine baseline differences. We used the t‐statistic to identify consistently altered loci between tumors and normal tissue. A list of the top 100 gene associated fragments that are significantly differentially methylated between tumors and normal tissue is provided in Supplemental Table 1 and was employed in the hierarchical clustering described in the following section. Some of our results confirmed previous findings for specific genes, such as the CpG island proximal to RUNX3 (p‐value = 8.4 e−7) and PITX2 (p‐value = 2.4 e−29) which have been shown to be inactivated in breast cancer (Maier et al., 2007; Harbeck et al., 2008; Subramaniam et al., 2009). We also identified several strong novel candidate biomarkers for breast cancer. Additionally, we found evidence that methylation of CpG islands upstream of certain microRNAs is correlated with tumorigenesis in this sample set. The CpG islands within 5 kb of miRNAs miR196a1 (p‐value = 0.00058), miR335 (p‐value = 2.2e‐14), miR124a3 (p‐value = 1.4e‐9) and miR423 (p‐value = 4.19e‐13) all show increased methylation in tumors when compared to normal tissue. Expression of one of these, miR335, has been previously shown to be lost in a majority of primary breast tumors and this loss of expression leads to an increased likelihood of metastasis (Tavazoie et al., 2008). The methylation status of three such candidates (GPR10, DRD5, CDKN1C) in a set of normal tissues and breast tumors is plotted in Figure 1a. These three candidates were then evaluated for significant methylation in tumor tissue using an independent sample set using bisulfite sequencing. All three candidates showed significantly increased methylation in tumor tissue while maintaining minimal methylation in matched normal tissues ( Figure 1b).

Open in a separate window

Figure 1

(a) Differentially methylated loci between tumors and normal tissue of 3 loci proximal to genes GPR10, DRD5 and CDKN1C. The samples are grouped by expression subtype (b) Bisulfite sequencing of selected loci in independent sample set validates significant methylation in tumors compared to normal tissue.

We next investigated DNA methylation in the clinical context of breast cancer subtypes, histology and prognosis. The 108 tumors with known expression‐based subtypes were divided into a discovery set (83 samples) for clustering analysis and a validation set (25 samples). We determined the features to be used for clustering in two steps consisting of non‐overlapping sets. The top 500 loci that varied most by standard deviation across 83 tumor samples were chosen to maximize the methylation diversity in the set (Supplemental Table 2 ). Conversely, the 100 loci described above (Supplemental Table 1) that distinguished tumors from normal tissue contain little information concerning subtype diversity, but serve as unifying characteristics that are common to all breast cancers.

We used these 600 features to cluster the 83 tumors and 11 normal samples for which the expression subtype data was available. Hierarchical clustering of the samples based on these six hundred loci gave us clustering that is remarkably similar to the one produced by expression analysis. Figure 2a shows the clustering of samples based on methylation and overlays the known expression cluster of each sample. As expected, the normal breast tissue samples clustered tightly. The expression subtype was determined by identifying which of the five centroids each sample correlates most to as described by Sorlie et al. (Sorlie et al., 2003). Cluster I was dominated by samples with a high correlation to the Luminal A centroid and anti‐correlation to the Basal‐like centroid. This is in contrast to cluster II where the samples were highly correlated to the Basal‐like centroid and anti‐correlated to the Luminal A centroid. Samples in cluster III were dominated by low correlation to each of the centroids. Most of the Luminal A samples were in cluster I (22/30 samples) and the Basal‐like samples were in cluster II (7/8 samples), cluster III had a mixture of subtypes, but 5 of 10 Normal‐like samples were in this low correlation class. We computed the absolute of the difference between the correlation values to the Luminal A centroid and the Basal‐like centroid as a measure of confidence in assignment of the expression subtypes. This confidence measure is significantly greater in the samples which cluster together in the methylation‐based clustering (Methylation clusters I and II) than in samples in the third methylation cluster using a Mann–Whitney test (p‐value = 0.0013). Finally, changing the number of loci used from top 500 most variant to top 250 or top 1000 did not significantly alter the clustering and retained the separation of the luminal A samples from basal‐like/ErB2+ samples (data not shown).

Open in a separate window

Figure 2

(a) Unsupervised clustering of breast tumors using methylation data groups samples similarly to those obtained by expression based analysis. Luminal A (Blue) subtypes (Cluster 1) form a distinctly separate group from the basal‐like(red) and ErbB2+(purple) subtypes. ER status, TP53 mutational status and Her status by FISH are plotted under the dendrogram (Black = positive, Gray = Negative, white = NA). Correlations to Luminal, Basal and ErbB2+ Centroid derived from expression data are plotted under the heatmap. Clusters I and II have unambiguous assignment of expression subtype, while cluster III contains samples with samples having weak correlations to multiple centroids. (b) fragment subset clusters luminal and basal groups. (c) Addition of 25 validation samples to the clustering retains the separation of luminal subtypes from basal/ErbB2 enriched subtypes.

Open in a separate window

Figure 2

(Continued)

We then investigated the possibility of improving the clustering based on expression subtypes using a selected subset of loci. We used a semi‐supervised version of our evolutionary search tool, a Genetic Algorithm (GA) based feature selection method (Schaffer et al., 2005). In this approach, we used the GA to select and evaluate subsets of loci for uniformity of the cluster members' expression subtype. This subset selection resulted in many smaller subsets outperforming the original 600‐loci set as a stratification tool. We show in Figure 2b one such subset which included just 45 of the original 600 loci. This subset substantially improves clustering of the subtypes, most notably the basal samples. Finally we added 25 new samples as validation into the clustering algorithm and we were able to cluster 9/9 Luminal A samples in the Luminal cluster, 9/10 Basal‐like/ErbB2+ in the Basal‐like/ErbB2+ cluster and 5/5 Normal‐like samples in the Normal‐like cluster (Figure 2c).

Interestingly, these chosen loci were not related to the genes that were identified in the expression subset as intrinsic gene set. Of the 500 most varying loci, 146 loci were found to be within 5000 bp of the transcriptional start site of a known gene. Many of these genes have been implicated in breast cancer, notably the homeobox gene clusters (HOXA2, HOXB13, DLX3), cell surface receptors (EGFR, the WNT family receptor, FZD8, toll‐like receptor, TLR2), keratins (KRT7) and forkhead proteins (FOXF1). Supplemental Figure 1 shows a subset of these genes and their functional categories.

The remainder were not near any annotated gene transcriptional start site. We looked at the evolutionary conservation and regulatory potential for these loci to evaluate the possibility that these could have any unannotated regulatory functions. Using data from the multi‐species conservation derived regulatory potential of the genome (Taylor et al., 2006), we compared the loci that possess discriminatory power to cluster breast cancer subtypes over a randomly chosen subset. These islands that contribute to clustering turned out to be not significantly more conserved than the rest of the genome ( Supplemental Figure 2 ). The absence of any clear increased conservation or regulatory regions in loci that have subtype‐specific methylation is perhaps surprising. The CpG islands which vary the most among tumors are clearly able to differentiate between Luminal and Basal‐like tumors. This observation leads us to the possibility that the etiology of the different tumor subtypes maybe reflected in the genome‐wide methylation patterns. When we investigated the possibility that the luminal and basal subtypes have different global methylation levels, we determined that the overall methylation state in each sample was comparable. We find no evidence of a “methylator phenotype” in breast cancer as has been seen in colorectal cancer (Toyota et al., 1999a; Toyota et al., 1999b). The overall levels of methylation (as determined by the number of loci in each of the “−1”,“0” and “+1” states) in luminal and basal tumors remain the same, but distinct patterns of methylation separate the two types.

We then specifically looked for genes that are differentially methylated in the various subtypes. In this analysis we used a t‐test to identify fragments that were most different between Luminal A expression subtype and the Basal‐like or ErbB2+ expression subtypes. We identified over 637 loci that were significant (p‐value<0.01 after correcting for multiple‐testing), including 360 loci that were mappable to a gene using our previously defined criteria (Supplemental Table 3 ). Interestingly, 56% of the significant loci were found to be near genes, compared to only 30% in the most varying loci used in the unsupervised clustering analysis. A histogram analysis of the distance of these fragments to the transcriptional start site of their associated genes showed that most of the significant fragments were within conventionally defined promoter regions with more than 50% of our significant genes within 500b of the TSS. ( Supplemental Figure 3 ). The table includes many gene families implicated in cancer and differentiation – HOX A family (HOXA2, HOXA6, HOXA7, HOXA10, HOXA11, HOXC10, HOXC5, HOXD13), FOX family (FOXC1, FOXD4, FOXD4L1, FOXD4L3, FOXF2, FOXP4, FOXQ1), growth factors and growth factor receptors (EGFR, FGF9, FGF19), matrix proteins (COL14A1, COL16A1, COL7A1, CLDN10, CLDN5) and protocadherins (PCDH10, PCDHAC2, PCDHGA10, PCDHGA11). The HOX genes have been previously described to have aberrant methylation in a variety of cancers, with the HOX A cluster aberrantly methylated in breast cancer (Novak et al., 2006). Indeed, we find many HOX genes to be significantly altered in their methylation status. Importantly when we looked at the subtypes of cancer, the HOXA gene alterations were more common in tumors belonging to the Luminal A expression subtype than the others. Increased methylation of HOXA2, HOXA7, HOXA10 and HOXA11 was observed in a majority of Luminal tumors, while their levels stayed similar to normal breast tissue in the other subtypes ( Figure 3a). A heatmap of the average methylation levels of a few interesting candidates from the list of differential methylated genes with subtype‐specific methylation are shown in Figure 3b. Interestingly, a majority of the significant loci have higher levels of methylation in luminal samples when compared to the basal‐like/ErbB2‐plus samples.

Open in a separate window

Figure 3

(a) Deregulation of methylation at HOX A genes in Luminal cancer subtypes (b) Heatmap of average methylation values of 100 selected genes groups by subtype. Transcription factors belonging to Homeobox family (HOXA, IRX, LHX1), forkhead family (FOXC1, FOXD4, FOXP2), cell adhesion molecules (KRT7, COL7A1, COL14A1, COL16A1) and protocadherins (PCDHAC2, PCDHGA10) have increased methylation levels in luminal tumors when compared to basal and ErbB2 subtypes.

2.1. DNA methylation loci as prognostic factors

One of the critical factors in managing treatment for breast cancer is identifying those at most risk for metastatic disease. We used distant disease free survival (DDFS) as the clinical end point to measure tumor aggressiveness. The classical clinical prognostic factors such as node status and tumor size significantly stratified patients into good and poor prognosis groups, while receptor status and adjuvant therapy status did not possess significant prognostic value. The distributions and prognostic value of the clinical variables are summarized in Supplemental Table 4. We then used this clinical data to identify genomic loci whose methylation status correlated with tumor aggressiveness. Loci were first filtered by choosing only those that were likely to be informative by requiring that the minority methylation state of a given locus contains at least 15 samples. This resulted in a total of 34,371 loci, of which 11,459 loci were found to be within 5000 bp from a transcriptional start site and 22,912 loci were not near any transcriptional start site.

We evaluated each of the above loci for their ability to stratify patients into good or poor prognosis groups depending on their methylation status. Using the Kaplan–Meier estimator of survival as the underlying statistical model, we identified the association of methylation states with DDFS. The Kaplan–Meier estimator calculates the probability of no systemic recurrence at a given time by using the time to systemic recurrence from retrospective data. Survival curves were estimated for each state of a given locus and the log‐rank (Mantel‐Haenzel) test was used to determine whether the two survival curves were significantly different from each other. We estimated statistical significance based on 1000 permutations of the time to distant metastasis data and chose loci that achieved a significance level less than 0.05 after multiple‐testing correction. This led to a total of 2559 loci chosen as significantly stratifying the patients into good and poor prognosis groups.

To find which of the above loci are providing prognostic information independent of other clinico‐pathological variables, we performed multivariate Cox regression analysis using other significant clinical variables. Of the 2559 significant loci, a total of 921 loci remained significant when included in a multivariate Cox regression with the other clinical variables. Of these 921 loci, 490 were found to be within 5,000bp from a transcriptional start site while 431 were deemed intergenic. It is notable that the ratio of genic to intergenic loci amongst the prognostic factors is significantly different compared to the ratio of the complete list of 34,371 loci included initially for survival analysis (p‐value < 10⁻³). Gene ontology enrichment analysis of the 490 genic loci revealed significant (p‐value<0.01 after multiple‐testing correction) enrichment of genes related to transcription factor activity, regulation of MAP kinase activity, cell proliferation, cell death, angiogenesis and neuronal development (Beissbarth and Speed, 2004). This suggests that loci associated with poor prognosis are more likely to be functionally important in the metastasis cascade, thus opening up the possibility of identifying potential biomarkers for prognosis as well as targets for treatment.

Kaplan–Meier Survival curves are plotted for a set of interesting candidates in Figure 4. These genes are involved in cancer‐related molecular functions such as cell death (TOP1 p‐value = 5e‐04; SEPT4 p‐value = 6e‐03), cell cycle (UBE2C p‐value = 9e‐04; MAPK12 p‐value = 1e‐03; TUBB3 p‐value = 3e‐02), cell fate commitment (KLF4 p‐value = 9e‐03; FGF12 p‐value = 4e‐03) and transcription factor activity (ELK1, p‐value = 5e‐03). The list of prognostic genes organized according to their significant gene ontology terms are presented in Table 2. Notably, the methylation state of a gene amongst the samples with poor prognosis is different from its state in normal breast tissue samples. The complete list of genes and intergenic loci whose methylation status correlated significantly with relapse likelihood is given in Supplemental Table 5.

Open in a separate window

Figure 4

Loci whose methylation status predicts likelihood of relapse.

Table 2

Hazard ratios of selected fragments along with their GO category.

Significant Gene Ontology Term (p‐value)	Gene Name	Methylation State in Normal Tissue	Poor Prognosis Methylation State	Significance of Survival Difference(p‐value)	Significance of Multivariate Cox Coefficient (p‐value)	Multivariate Hazard Ratio (lower 0.95, upper 0.95)
Cell Death(1e‐03)	TOP1	M	UM	2e‐05	5e‐04	3.677 (1.75, 7.72)
	TDGF1	UM	M	6e‐05	9e‐04	4.08 (1.77, 9.38)
	CIDEB	UM	M	5e‐05	2e‐03	4.02 (1.63, 9.89)
	UNC5A	UM	M	4e‐03	0.027	2.74 (1.12, 6.7)
Cell Fate Commitment(6e‐04)	SPRY1	M	UM	1e‐04	2e‐03	4.0 (1.65, 9.73)
	FGF12	UM	M	1e‐04	4e‐03	3.31 (1.47, 7.49)
	KLF4	M	UM	1e‐03	9e‐03	2.63 (1.27, 5.46)
	OLIG2	UM	M	5e‐03	0.023	2.56 (1.14, 5.76)
Cell Proliferation(6e‐08)	KI‐67	M	UM	5e‐03	0.004	3.43 (1.48, 7.93)
	HDAC4	M	UM	1e‐03	0.016	2.62 (1.19, 5.75)
Cell Cycle Process(7e‐04)	KIF2C	M	UM	2e‐03	8e‐03	2.65 (1.28,5.49)
	UBE2C	M	UM	1e‐05	9e‐04	3.57 (1.68, 7.59)
	MAPK12	M	UM	9e‐04	1e‐03	3.51 (1.60, 7.67)
	TUBB3	M	M	9e‐03	0.035	3.19 (1.09, 9.39)
Vasculature Development(8e‐04)	VEGF	M	UM	8e‐03	7e‐03	3.01 (1.34, 6.75)
	KLF5	M	UM	4e‐03	4e‐03	3.56 (1.49, 8.53)
	BTG1	M	UM	2e‐03	0.013	2.65 (1.22, 5.74)
Transcription Factor Activity(1e‐04)	LHX5	UM	M	6e‐03	0.032	2.94 (1.10, 7.87)
	LHX2	M	UM	4e‐03	0.028	2.32 (1.09, 4.91)
	ONECUT2	M	UM	8e‐03	0.029	2.74 (1.11, 6.79)
	ONECUT1	UM	M	8e‐03	0.023	2.8 (1.15, 6.83)
	FOXH1	UM	M	8e‐05	7e‐05	5.09 (2.28, 11.36)
	ELK1	UM	M	1e‐03	5e‐03	3.53 (1.46, 8.53)
	JUN	M	UM	4e‐03	0.025	2.44 (1.12, 5.33)
	GSC	M	UM	3e‐04	2e‐03	3.69 (1.61, 8.45)

Open in a separate window

We find demethylation of the Topoisomerase I (TOP I) promoter increased the likelihood of relapse in those patients. Other candidate genes include Goosecoid (GSC, p‐value = 3e‐04), which has been previously shown to promote metastasis through its role in epithelial‐mesenchymal‐transitions (Hartwell et al., 2006) and whose demethylation is found to be associated with poor prognosis in our study; demethylation of Vascular endothelial growth factor, VEGF (p‐value = 8e‐03), whose role in angiogenesis is well known is also associated with poor prognosis; and ONECUT1 (p‐value = 8e‐03), whose methylation has been previously associated with cervical cancer and is associated with poor prognosis in our study. We found that the demethylation of TP53BP2 (p‐value = 1e‐03), a protein that binds to TP53 leads to significantly poorer prognosis independent of treatment and other clinical variables. Due to its reported role in regulating the apoptotic function of TP53 (Sullivan and Lu, 2007), we investigated the relationship of TP53BP2 methylation and relapse risk in TP53 wild type and TP53 mutated populations. Interestingly the prognostic value of TP53BP2 methylation status was completely eliminated in TP53 mutated populations, revealing an effect modifying link between TP53 mutation and TP53BP2 methylation based prognosis ( Figure 5 ).

Open in a separate window

Figure 5

(a) Methylation status of TP53BP2 is associated relapse risk, (b) risk associated with TP53 mutation status (c) BP2 methylation status is informative of relapse risk only in patients with wild type TP53 but not (d) TP53.

2.2. Effects of methylation on Gene expression

The OMS study samples have been profiled using cDNA arrays to measure RNA levels in each sample relative to a pooled RNA reference and expression data are publicly available. The absence of absolute expression levels precluded a sample by sample correlation to the absolute methylation levels identified in this study. However, it is possible to pool the samples according to their methylation states and identify the difference in relative expression levels between these groups. For each locus, we grouped samples according to their methylation state (−1,0,+1). We then tested if there is a difference between expression levels of a given gene between 2 methylation states using the one‐sided Wilcoxon signed rank test. We used a definition of a fragment being within 5 kbp upstream to 2 kbp downstream of a transcriptional start site to identify that fragment as being associated with the gene. 9487 unique genes from the expression data set had at least one fragment mapped to this region. We analyzed the effects of methylation of these loci on expression levels of associated genes. Of these, 8393 genes had at least 10 samples exhibiting at least 2 methylation states thereby enabling us to compare relative expression level differences with methylation state. 2853 (33%) genes showed a significant anti‐correlation of the expression levels to the methylation state of the corresponding fragment (p‐value<0.05). 146 of the 500 fragments used in the unsupervised clustering were gene associated by criteria described above. Of these 146, 79 were present in the expression data. 33 of those 79 genes showed significant anti‐correlation of gene expression with methylation. Additionally, 313 genes that were significant in the survival analyses had expression data and 137 of those were found to have significant anti‐correlation between gene expression and methylation status. Of the 362 genes with significant differential methylation between basal‐like/ErbB2+ and luminal subtypes, 200 could be analyzed for correlation to expression. 118 of these genes showed significant anti‐correlation of gene expression levels to methylation states.The mean expression values of the samples in each methylation state are plotted for 50 genes in a heatmap in Figure 6.

Open in a separate window

Figure 6

Correlation of expression levels to methylation states of the associated fragments. Heatmap of the mean of the expression values for a selected list of 50 genes with significant anti‐correlation to methylation states is plotted. White colors indicate no samples were present in that methylation state for the gene. The number of individual samples in each state varies by gene and thus is not plotted.

3. Discussion

In this study, we performed genome‐wide scans of CpG island methylation patterns in over 100 breast tumors and normal breast tissues. The primary finding of our analysis is that luminal breast tumors have different methylation profiles when compared to other breast tumor subtypes. The pattern of methylation in Luminal tumors is distinctly different from tumors of basal‐like or ErbB2+ origin. We find that this difference is reflected throughout the CpG islands in the genome and is not limited to functional genes. We compared our results with those of the parallel study of Rønneberg et al using the Illumina Golden Gate array and reported concurrently along with this study. Our own clustering results produced three clusters (Figure 2c). Cluster I consists of mostly Luminal A subtypes, Cluster II has ErbB2+ and basal subtypes and cluster III has normal samples together with Normal‐Like and some Luminal subtypes. We determined that 72 samples were analyzed in both studies. The Rønneberg et al study also identified three methylation based clusters. Two of these clusters consisted of a majority of Luminal A expression subtypes (Clusters I and III) and one with Basal‐like and ErbB2+ subtypes (Cluster II). There were significant additional differences in ER and TP53 mutational status between Cluster I and II but not with the Cluster III. When we looked at the overlaps between our Luminal A cluster (Cluster I) with those of Rønneberg et al., 32/34 samples were identified in the Luminal majority clusters (Clusters I and III in Rønneberg et al.). Additionally, 23/28 samples in Cluster II in our study were also identified as belonging to cluster II in the Rønneberg study. These overlaps provide strong evidence of the distinct methylation differences between Luminal A subtype and the basal/ErbB2+ subtypes. Our results predict that there are fundamental changes occurring in luminal tumors that are then reflected in the genome‐wide methylation patterns. Interestingly, luminal origin tumors had significantly more methylation in fragments associated with transcription factors implicated in development and differentiation – notably the HOX A and homeobox genes and forkhead family.

These findings complement recent reports on DNA methylation patterns seen in different cell types in breast tissue (Bloushtain‐Qimron et al., 2008). Bloushtain‐Qimron et al. identified methylation patterns from FACS sorted cells from normal breast tissue; CD24+ (Luminal), MUC1+ (Luminal progenitor), CD10+ (myoepithelial progenitor) and CD44+ (multipotent progenitor). They determined that CD24 + luminal cells have increased methylation at a number of development related transcription factors when compared to CD44 + basal progenitor cells. Our own data show significant overlaps of the methylation patterns in luminal subtypes with those found in CD24 + breast cells and patterns in basal subtype overlap with multipotent CD44 + progenitor cells. A comparison of differentially methylated genes between the Bloushtain‐Qimron study and our own results showed remarkable agreement – 10 of the top genes identified in their study are consistent with our results. Table 3 shows overlaps between our study and those of Bloushtain‐Qimron et al. These findings suggest that the genome‐wide methylation pattern of the tumors reflects the methylation pattern of the cell of origin. Supporting this hypothesis is the finding that FOXC1 and HOXA10, two of the four markers identified by Bloushtain‐Qimron et al. are among the top most differentially methylated markers in our study. During the preparation of this manuscript, two new studies, (Flanagan et al; Holm et al.) reported their findings on DNA methylation patterns in breast cancer. The study by Holm et. al. showed significant differences in the DNA methylation patterns of luminal A, luminal B and basal tumors using the Illumina Beadstudio Methylation Module. These results complement our own findings, although their platform covered only 807 cancer‐related genes compared to the genome‐wide CpG islands platform used in our study. Due to this limitation, they did not identify the novel findings in our own study that DNA methylation patterns in the breast cancer subtypes are not limited only to cancer‐related genes but reflected genomewide in both gene associated and non‐gene associated CpG islands. The study by Flanagan et al., determined DNA methylation profiles in 33 familial breast cancer samples. This study primarily focused on identifying differential methylation patterns driven by BRCA1 mutation status and did not report methylation patterns being associated with subtypes.

Finally we discovered candidate loci for cancer prognosis with molecular functions associated with cancer ‐ Cell Cycle and Proliferation (Ki‐67, UBE2C, KIF2C, HDAC4), vasculature development and angiogenesis (VEGF, BTG1, KLF5), transcription factors with roles in cell fate commitment (SPRY1, OLIG2, LHX2 and LHX5). Interestingly over‐expression of proliferation markers such as Ki‐67 and UBE2C are already being considered as prognostic markers (de Azambuja et al., 2007; Loussouarn et al., 2009) and correspondingly their methylation levels are lower in patients with poor prognosis. Topoisomerase I, which has previously implicated in chemotherapy resistance, was found to be also a risk factor for recurrence. The methylation of TOP1 may also have clinical implications for patients. 30% of breast tumors in our study have an unmethylated promoter and these patients were observed to have a higher likelihood of relapse Topoisomerase expression has been linked to the development of resistance to topoisomerase inhibitors (Burgess et al., 2008). Demethylation of TP53BP2 leads to significantly poorer prognosis (p‐value = 1e‐03) independent of treatment and other clinical variables. TP53BP2, which interacts with p53 has been implicated in gastric cancer (Ju et al., 2005). The existence of this effect in only p53 wild type tumors is intriguing and further study is needed to determine the functional nature of this link and likely mechanism of the interaction. An exhaustive analysis of all significant genes for interaction between themselves and p53 mutation status would provide further insight into the mechanisms of disease progression. Similarly, the interaction between all individual methylation loci which have prognostic ability have to be evaluated in future studies to identify combination markers which have better performance. The interaction between the individual loci may also provide insight into the different molecular pathways that are affected by treatment.

In summary, our study has found evidence for a strong association of DNA methylation within breast cancer subtypes. The question of whether this association is causal should be the focus of future studies – especially the high number of significantly methylated genes in luminal A subtypes in comparison with the basal or ErbB2+ subtypes. This is especially intriguing since we could not find any difference in the overall methylation levels among the subtypes. While this could be due to innate differences between cell types that may give rise to Luminal A and basal subtypes, there is no convincing evidence that points to the differences in the cell type of origin in the subtypes. A focused study addressing these questions would greatly advance our understanding of the breast cancer subtypes and their implications for treating patients with efficacy.

4. Methods

Supporting information

Supplementary material 1. 

Notes

Kamalakaran Sitharthan, Varadan Vinay, Giercksky Russnes Hege E., Levy Dan, Kendall Jude, Janevski Angel, Riggs Michael, Banerjee Nilanjana, Synnestvedt Marit, Schlichting Ellen, Kåresen Rolf, Shama Prasada K., Rotti Harish, Rao Ramachandra, Rao Laxmi, Eric Tang Man-Hung, Satyamoorthy K., Lucito Robert, Wigler Michael, Dimitrova Nevenka, Naume Bjorn, Borresen-Dale Anne-Lise, Hicks James B., (2011), DNA methylation patterns in luminal breast cancers differ from non‐luminal subtypes and can identify relapse risk independent of other clinical variables, Molecular Oncology, 5, 10.1016/j.molonc.2010.11.002. [Europe PMC free article] [Abstract] [Google Scholar]

Contributor Information

Sitharthan Kamalakaran, Email: moc.spilihp@khtrahtis.

Vinay Varadan, Email: [email protected].

Hege E. Giercksky Russnes, Email: ude.dravrah.icfd@senssur_egeh.

Dan Levy, Email: ude.lhsc@yvel.

Jude Kendall, Email: ude.lhsc@lladnek.

Angel Janevski, Email: [email protected].

Michael Riggs, Email: ude.lhsc@sggir.

Nilanjana Banerjee, Email: [email protected].

Ellen Schlichting, Email: [email protected].

Rolf Kåresen, Email: [email protected].

Man-Hung Eric Tang, Email: ude.lhsc@gnatm.

Robert Lucito, Email: ude.lhsc@oticul.

Michael Wigler, Email: ude.lhsc@relgiw.

Nevenka Dimitrova, Email: [email protected].

Bjorn Naume, Email: [email protected].

Anne-Lise Borresen-Dale, Email: on.oiu.muidar@bla.

James B. Hicks, Email: ude.lhsc@skcih.

References