Identification of Breast Cancer DNA Methylation Markers Optimized for Fine-Needle Aspiration Samples

This article is highlighted in the In This Issue feature, p. 2143

Random periareolar fine-needle aspiration (RP-FNA) samples are frequently used in phase 2 trials of breast cancer prevention for assessments of surrogate endpoint biomarkers (1–7). DNA methylation is a potentially reversible early event in breast carcinogenesis (8–11), and is readily assessable in paucicellular RP-FNA samples (12–15). For some genes, higher levels of methylation are identified in benign breast cells from cancer patients than similar samples from unaffected women (16–18) and there is, often, a correlation between specific methylation changes observed in benign breast tissue and the associated cancers (16, 19, 20). This suggests that measures of DNA methylation in benign breast tissue may predict breast cancer risk, and, in fact, DNA methylation of some genes is observed more frequently in breast epithelial cells obtained by RP-FNA from high-risk women as compared with lower risk women (17, 18, 21); in addition, tamoxifen has been shown to reduce methylation measured in RP-FNA samples for some genes (7). Although the risk signals generated by these candidate gene studies are intriguing, the clinical relevance of hazard ratios in the 2.5 to 3.0 range is marginal at best.

Clinically useful RP-FNA-based breast cancer risk stratification may be achievable with better DNA methylation markers. Methylation arrays have greatly expanded the list of potential markers (22–28), but have not necessarily identified genes whose expression is regulated by promoter region methylation. We used a genome-wide 5-aza-2′-deoxycytidine (5AZA) approach to identify genes whose expression is upregulated with demethylation. The gene lists were subsequently refined based on several phases of qualitative and quantitative methylation-specific PCR (MSP) assays targeting specific promoter regions in benign and malignant breast samples and lymphocytes. A selected subset of markers was assessed in a large panel of clinical breast cancer FNA and benign breast RP-FNA samples to confirm clinical relevance.

Six breast cancer cell lines were selected for demethylation studies based on known tumor suppressor gene expression regulation by promoter region hypermethylation: HCC1569 (CCND2), HCC1954 (SCGB3A1, APC, RASSF1A), MCF-7 (RAR-β2), MDA-MB-231 (ESR1), UACC3199 (BRCA1), and BT-549 (hypermethylator phenotype). Apart from MCF10A, we specifically avoided immortalized benign human mammary epithelial cell lines for this experiment as these cells frequently show tumor suppressor gene methylation (e.g., p16) and gene expression profiles that are intermediate between normal breast epithelial cells and breast cancer (29). Instead, we opted to test six first-passage benign human mammary epithelial cell cultures (HME) generated in serum-free media from small fragments of normal breast tissue obtained from young women undergoing fibroadenoma excision. Fresh benign breast tissue was minced and digested overnight on a rotator at 37°C in 1 mg/mL collagenase (C2674, Sigma-Aldrich) in Mammary Epithelial Cell Growth Medium (MEGM; NC9523177, Fisher Scientific) and then centrifuged at 200g for 3 minutes. The pellet was rinsed with PBS, centrifuged again, and then treated with 0.25% trypsin in EDTA at 37°C for 3 minutes to dissociate the cell clumps. After neutralizing the trypsin, the cells were centrifuged for 3 minutes at 200g, resuspended in MEGM with 5% FBS, and plated into tissue culture flasks. The media was changed to MEGM without serum the next day. Contaminating fibroblasts were removed by a quick trypsinization 3 to 5 days later.

Dose-dependent cytostasis was observed in 5AZA-treated cells, and was most pronounced for the primary benign breast epithelial cell cultures. The 5AZA dose (0.5 μmol/L) was selected based on evaluation of growth curves and induction of BNC1, SERPINB, and TKTL1 gene expression measured by reverse transcription polymerase chain reaction (RT-PCR) in benign and malignant cells. The breast cancer cell lines, HME cultures, and MC10A cells were treated with 0.5 μmol/L 5AZA (Sigma-Aldrich) in dimethylsulfoxide (DMSO) or DMSO alone for 6 days, after which RNA was prepared using the Illumina TotalPrep kit (AMIL1791, Life Technologies). This RNA was used for whole-genome expression analysis.

Whole-genome expression was assessed using the Illumina HumanWG-6-v3 chip. The initial analysis was done in GeneSpring as follows. Flags were set to present for entities with detection P values of 0.2 or less and absent for entities with detection P values of more than 0.4. Raw signals with values less than 1 were reset to 1. Quantile normalization was used, and baseline transformation performed on the basis of the median expression level for all entities. Of the 48,701 entities included on the Illumina array, flags were present or marginal for 75% of the samples for 27,324 entities, and these were retained. The complete dataset is available on the NCBI GEO site (GSE41692).

Publically available informatics tools were used to limit the gene list identified in the 5AZA experiment to those most likely to yield useful epigenetic markers (30–34). Specific exclusions included duplicate genes (192), no CpG island (87), known imprinted genes (5), X-chromosome genes and likely imprinted (30), pseudogenes or noncoding RNA (6), poorly annotated or discontinued records (12), and CpG island more than 2,000 bp proximal to the transcription start site (33).

Qualitative methylation-specific PCR (MSP) assays were designed and run as described previously (21). For genes previously shown to be epigenetically regulated, published sequences were used. For novel genes, primers were designed within the CpG island as close to the transcriptional start site as possible (median 5′ distance from the transcription start site was 97 bp with a median product size of 113 bp). Primer sequences used in the initial MSP screen are shown in Supplementary Table S1. For the second MSP screen in clinical samples, HotStar taq (Qiagen) was used, and MgCl₂ concentrations and PCR cycles optimized to permit the detection of 60 pg of methylated NaBS-treated gDNA in a 6-ng DNA sample (Supplementary Table S2).

Quantitative multiplex MSP (QM-MSP) assays were designed as described previously by Fackler (35) using the same conditions and quality assurance standards we have previously described (17, 18). Results were expressed as a methylation fraction, which is calculated as methylated copies/(methylated copies + unmethylated copies). Primer and probe sequences are shown in Supplementary Table S3.

Five markers were selected from the QM-MSP screen for assessment in an expanded clinical sample set that included all of the samples used in the initial QM-MSP screen with the addition of new, prospectively acquired samples. This research was performed in accordance with an assurance filed with and approved by the U.S. Department of Health and Human Services. Institutional Review Board approval was obtained, and informed consent was documented in writing for each participant. Unselected patients with untreated stages 1 to 3 primary breast cancer (n = 52) underwent FNA sampling of the primary tumor at the time of definitive breast cancer surgery. In addition, benign breast epithelium was obtained by RP-FNA from patients with breast cancer (n = 52) and women never diagnosed with breast cancer (n = 90) as previously described (18, 21). Samples not meeting quality assurance standards for QM-MSP (∼15%) were excluded from analysis; therefore, the number of samples shown in the tables is lower than the number collected. The total number of samples included in this analysis (archival plus prospectively acquired) is 97 breast cancers (one bilateral breast cancer), 104 benign RP-FNA samples from cancer patients, and 223 benign RP-FNA samples from women never diagnosed with breast cancer.

Clinical information was prospectively collected on case report forms and tumor data, including tumor size, nodal status, associated ductal carcinoma in situ (DCIS), expression of estrogen receptor (ER), progesterone receptor, and Her-2/neu, proliferative index (Ki67), and p53 expression, were subsequently abstracted from the pathology reports.

For the purposes of summarization and comparison, quantitative methylation data are typically dichotomized using an arbitrary threshold value above which a sample is scored as methylated and below which a sample is scored as not methylated. This approach artificially constrains the data, reducing it to the level of qualitative MSP data. To avoid this when making comparisons between patient groups, we dichotomized the patient variable of interest, ER status for example, and then plotted the proportion of variable (+) and variable (−) cases that exceeded a continuous range of possible methylation threshold values (e.g., 10⁻⁸ to 0.5). These two curves (Fig. 1) were compared as follows. Methylation values were randomly permuted to destroy the association, if any, between methylation of a specific gene and the selected clinicopathologic variable of interest. For each permuted data set, we calculated the proportions of cases in each of the two categories as in Fig. 1. Note that if there is no association between the methylation and the clinicopathologic variable, the two curves are expected to be identical. We then computed the differences between the two curves and used their summation as a measure of the similarity between the two groups. A null distribution for the sum of the differences is developed over the course of 10,000 iterations, and P values are calculated based on the null distribution. To account for the multiple comparisons included in this analysis, false-discovery rate (FDR)–corrected P values were also calculated by the method of Benjamini and Hochberg (36).

Initially, 5AZA-induced gene expression was used to identify potential breast cancer-relevant methylation markers. The resulting gene list was refined using informatics filters followed by three phases of MSP testing in panels of cell lines, benign primary breast epithelial cell cultures, lymphocytes, and clinical samples. A subset of five markers fully optimized for RP-FNA samples was tested in an extended clinical sample set to confirm clinical relevance. The marker discovery and validation pipeline is shown in Supplementary Fig. S1.

Volcano plots initially identified 72 genes whose expression was induced 2-fold or more with an FDR of less than 0.05 in breast cancer cell lines treated with 5AZA. Expression of two genes was reduced by 5AZA by 2-fold or more with FDR less than 0.05 (HDAC4 and SMYD3). Hierarchical cluster analysis of baseline expression of 72 genes induced by 5AZA discriminated perfectly between the benign and malignant samples (Fig. 2). The pattern of expression for these 72 genes was highly consistent among the benign breast epithelial cells and more variable among the cancers. Although MCF10A cells clustered with the benign breast epithelial cells, they were recognized as a distinct subclass. Thirty-six of these 72 genes showed reduced baseline expression in the breast cancers as compared with the benign cells (Fig. 2, orange box), but only 20 of these passed the informatics filters required for advancement to the MSP screens. To reduce the probability of overlooking useful genes, we ultimately advanced genes with 4-fold or greater induction by 5AZA in one or more breast cancer cell lines, or 1.5-fold or greater induction by 5AZA in two or more breast cancer cell lines. In total, 645 entities, representing 453 unique genes, were selected from the 5AZA experiment. Of these, 280 passed the informatics filters, but MSP assays could be optimized for only 271 of these.

There were 31 (11.4%) transcription factors represented among the 271 5AZA-induced genes, which is a significant enrichment over the 1,857 (7.4%) transcription factors included among the 25,186 unique Entrez IDs found in the Illumina genome (P = 0.011). This list included only three homeobox genes, HOXB5, HOXD1, and MSX1, which is consistent with the total number of homeobox genes found in the Illumina genome (P = 0.974). Gene ontology analysis showed enrichment for several biologic processes, including response to stimulus (corrected P = 0.013), multiorganism processes (corrected P = 0.006), and developmental processes (corrected P = 0.006). Enrichment for response to stimulus genes was most significant for the response to stress terms (corrected P = 0.0008).

The 271 genes identified in the 5AZA experiment were tested by qualitative MSP assays in 10 human breast cancer cell lines, six benign early passage primary HME cultures, and four lymphocyte samples (Supplementary Table S1). These breast cancer cell lines were selected to include ER-positive, ER-negative, HER2-positive, and HER2-negative cancers and did not overlap with the cell lines used in the 5AZA experiment. Because our final implementation is in RP-FNA samples, many markers underwent more extensive testing in lymphocytes. Markers generating any signal in lymphocytes were excluded prior to validation in clinical samples. Ninety-four of these 271 genes were methylated in breast cancer cell lines but not in benign HME cultures (34 of these genes were also methylated in lymphocytes). Of these, 72 did not occur on previously published methylation array-based marker lists (Fig. 3; refs. 22–25, 37). Many common breast cancer methylation markers were absent from our list, including 13 that we have recently been assessing for breast cancer–risk stratification: LOX, BNC1, CFTR, HS3ST2, CTSZ, ESR1, BRCA1, CDH13, RASSF1, CCND2, APC, RAR-β2, and SCGB3A1. Nine of these 13 genes did not meet expression array detection P value thresholds, and the others failed to meet thresholds for 5AZA-induced upregulation. Notably, only six of these 13 genes were among the 1,136 markers identified by Van der Auwera using the Infinium Human Methylation-27 platform (22), and only two were among the 263 identified by Hill using the same platform (24). None occurred on any of the other array-derived lists shown in Fig. 3.

From the 271 5AZA-induced genes assessed by qualitative MSP in breast cancer cell lines and HME cultures, 102 were excluded because of methylation signals in lymphocytes, 99 because no methylation signals could be detected for any cell type, and 21 because of very infrequent methylation in the cancer cell lines. The remaining 49 markers were assessed by qualitative MSP in a panel of fresh-frozen primary breast cancers (15), benign RP-FNA samples from untreated women recently diagnosed with breast cancer (5), and benign RP-FNA samples from women never diagnosed with breast cancer (10). Methylation prevalence was 20% or greater for cancers than benign samples for 33 of these markers and 40% or greater for 15 markers (Supplementary Table S2). Eight of these later markers do not occur on the methylation array-based marker lists included in Fig. 3 (HLA-B, LAT2, FBLN2, VCAN, ADM, FLNC, ARTN, and CLDN1; refs. 22–25, 37).

The TCGA data presents an opportunity to understand the characteristic of our 5AZA-derived marker list in relation to methylation array-derived lists. TCGA data is derived from the Infinium Human Methylation-450 array, which quantifies the level of methylation (β-value) for each of 485,000 CpGs distributed across the genome for 627 breast cancers and 97 benign breast samples. We initially identified 140 loci associated with unique named genes where the mean β was 0.4 or higher for malignant than benign samples (Fig. 3). Of the 33 5AZA genes differentially methylated in cancer as compared with benign samples (tumor minus benign methylation prevalence, ≥0.2), only one occurred among these 140 most differentially methylated TCGA genes. This gene was GNE. Next we expanded the TCGA list to include any probe with a mean β difference greater than 0 for a unique named gene (n = 9,509). This list was ordered by FDR for tumor versus benign differentiation, to determine how many of our 33 differentially methylated genes could have been derived from TCGA data and how far down the list they occur. Fifteen of our 33 differentially methylated 5AZA genes were recognized as differentially methylated in the TCGA data, and their positions on the sorted list ranged from 199 for AKR1B1 (FDR = 8.5 × 10⁻¹³⁵) to 6721 for GPX1 (FDR = 1.3 × 10⁻¹⁹). Finally, we identified the TCGA probe closest to our amplicon for the 49 genes assessed by MSP in clinical samples and compared the tumor versus benign discrimination from the MSP data (difference in methylation prevalence) with the TCGA data (difference in mean beta). We were able to identify TCGA probes 2 to 279 base pairs (mean, 71) from the center of our MSP amplicons for each of these 49 genes. There was significant correlation between tumor versus benign discrimination for MSP and TCGA for these 49 markers (Spearman r = 0.406, P = 0.004). In addition, the difference in mean β values for malignant versus benign samples from the TCGA data was 0.085 for the 15 markers with the greatest difference in methylation prevalence by MSP (≥0.4 difference threshold) as compared with only 0.028 for the 16 least discriminatory markers (<0.2 difference threshold, P = 0.025). This suggests that the best 5AZA-derived markers would be identifiable within the TCGA data if the filtering algorithms were generously relaxed.

A subset of the genes with higher methylation prevalence in malignant than benign clinical samples was selected for QM-MSP testing. In order to select genes whose methylation status was largely independent of the methylation status of other genes, hierarchical clustering of MSP methylation results for 15 breast cancers was performed (data not shown). This identified two gene clusters. The first cluster included genes that were methylated in the majority of the cancers (AKRB1, PSAT1, CYP24A1, GBP4, HLA-B, LAT2, CPNE8, and UCHL1) and the second cluster included genes methylated in only a fraction of the cancers and clustering with one to five other genes. Four markers from clusters 1 and 8 distributed across cluster 2 (12 markers total) were selected for QM-MSP validation in a panel of 51 archival primary breast cancer FNAs, 59 archival benign RP-FNA samples from untreated women recently diagnosed with breast cancer, and 145 archival benign RP-FNA samples from women never diagnosed with breast cancer. These genes included WDR66, UCHL1, FBLN2, PSAT1, CLDN1, BIRC3, CYP24A1, CCNA, GNE, CXCL14, CPNE8, and AKR1B1. Of these, all except BIRC3, GNE, and WDR66 were differentially methylated between malignant and paired benign samples with P value between 0 and 0.013 by sum-of-difference analysis (Supplementary Table S3).

Methylation prevalence for these 12 genes was compared between the 15 fresh-frozen primary breast cancers used in the initial MSP screen and the 51 breast cancer FNA samples assessed in the QM-MSP screen using greater than 0.1% as the threshold for classifying a sample as positive, as this is the reported sensitivity of MSP (38). Methylation was identified at a modestly higher frequency in fresh-frozen breast cancer tissue than in breast cancer FNA samples for each of the 12 genes (mean 0.53 vs. 0.43, P = 0.13), and correlation between the samples for methylation prevalence was poor (Spearman r = 0.343, P = 0.298). This may reflect tumor heterogeneity, which is missed by FNA sampling.

Four markers with higher levels of methylation in malignant than in paired benign samples were selected for QM-MSP validation in an expanded clinical sample set (CPNE8, PSAT1, CXCL14, and CLDN1). GNE was included as the fifth marker because it had shown a trend for differential methylation in malignant as compared with paired benign samples (P = 0.085), and methylation in benign samples had only been observed for women recently diagnosed with a primary breast cancer. The sample set included 97 breast cancer FNAs (one bilateral breast cancer, Table 1), 104 benign RP-FNA samples from patients with breast cancer, and 223 benign RP-FNA samples from women never diagnosed with breast cancer.

Each of these five markers is methylated at a greater frequency in primary breast cancer FNAs than in benign RP-FNAs, and each, except GNE, shows significantly greater methylation in primary cancers than paired benign samples (Table 2). The sum-of-differences permutation method was used to assess the association between methylation of specific genes and various clinical and pathologic features of the patients and tumors (Table 3). Most notably, the PSAT1 methylation was associated with low-grade, low-proliferation, hormone receptor (HR)-positive, lymph node–positive breast cancer in postmenopausal Caucasian women and with infiltrating lobular carcinoma. Conversely, GNE methylation was associated with high-grade, HR-negative breast cancer in younger women. Most of the markers showed greater methylation in tumors from Caucasian women than African-American women, and CXCL14 methylation was associated with HER-2/neu-positive breast cancer. Methylation of CLDN1 in a benign RP-FNA sample was marginally predictive of methylation of CLDN1 in the paired cancer sample (P = 0.058). This was not observed for the other four markers.

Whole-genome array-based approaches have identified hundreds of potential methylation markers in breast cancer (22–25, 37) using probes specific for individual CpGs found near named genes. Studies using the same methylation platform will generally identify different markers depending on the samples tested and the filtering criteria used. These marker lists will be even more divergent for different platforms (Fig. 3). Methylation arrays examine individual CpGs. One highly methylated CpG does not necessarily equate to dense CpG methylation in a region critical for regulation of gene expression. Indeed, for a subset of genes, CpG methylation is a consequence of transcriptional repression and not a cause as evidenced by failure of 5AZA treatment to induce expression (39).

In order to identify DNA methylation markers with a high probability of regulating gene expression, we started by identifying genes expressed at higher levels in benign breast epithelial cells than breast cancer cells and then selected genes that were induced by 5AZA treatment in the cancer cells. Two phases of MSP screening were used to identify 33 markers differentially methylated in clinical breast cancer samples as compared with clinical benign samples without generating signals in lymphocytes (Supplementary Table S2, δCA-B9 ≥ 0.2). These are breast cancer-relevant methylation markers optimized for FNA samples.

The clinical relevance of a subset of these markers was demonstrated in an expanded clinical sample set using QM-MSP. Consistent with other studies (19, 23, 24, 26, 28, 40–44), DNA methylation was observed more frequently in HR-positive than HR-negative breast cancer, and this difference was most marked for PSAT1, CPNE8, UCHL1, and AKR1B1. Notably, PSAT1 methylation was associated with low-grade, HR-positive, lymph node–positive breast cancer in postmenopausal women and infiltrating lobular cancer. In addition, GNE methylation was associated with high-grade, HR-negative breast cancer and CXCL14 methylation with HER2-positive breast cancer.

Only 12 of the 33 5AZA-selected genes we identified as differentially methylated in malignant as compared with benign clinical samples were also identified by methylation array studies (CYP24A1, NBL1, CPNE8, AKR1B1, UCHL1, HBA2, GSTP1, GBP4, PSAT1, IRF7, PYCARD, and GNE; refs. 22–25, 37). In addition, many common tumor suppressor genes well-documented to be regulated by promoter-region hypermethylation were missed by our 5AZA approach and by methylation array approaches (e.g., BRCA1, CFTR, CTSZ, ESR1, LOX, RAR-β, and SCGB3A1) although CCND2, APC, RASSF1, CDH13, BNC1, and HS3ST2 do occur among the 1,136 genes identified by Van der Auwera using the Infinium Human Methylation-27 platform. These genes were missed by our 5AZA approach because of failure to meet expression level or quality thresholds on the initial gene expression array as described in the Results section. Methylation array data are generally filtered based on mean β-scores for individual CpGs in order to retain genes where neighboring CpGs are methylated for the majority of DNA copies in the majority of samples. This will overlook genes that are methylated in a relatively small subpopulation of cells or methylated in only a small subset of tumors. Methylation-specific PCR, used for methylation marker discovery in the past and for initial validation in our study, will generate a signal when only 0.1% of DNA copies are methylated (38). This may permit recognition of minor cell populations important for tumor biology and maintenance. In this regard, the sum-of-difference permutation method we described for analysis of quantitative methylation data essentially compares the “sizes” of methylated populations between two groups and is ideal for recognizing differences in the relative abundance of minor cell populations.

The clinical utility of methylation markers in breast cancer has not yet been well established, and it is unclear whether the greatest value will come from loci methylated in major cell populations in most cancers, or from genes regulated by promoter methylation in minor cell populations. Our primary interest is in RP-FNA markers for breast cancer risk stratification. For this application, identification of minor cell populations (e.g., tumor-initiating cells) may be particularly relevant. However, our data suggest that tumor FNA samples may not capture the entire methylation spectrum identifiable in fresh-frozen tumor tissue. Additional work is required to recognize the markers most relevant to tumor initiation, to validate tumor FNA as an adequate sampling approach, and to map relevant methylation changes in benign breast tissue.

No potential conflicts of interest were disclosed.

Conception and design: D. Bu, C. Lewis, A.F. Gazdar, D.M. Euhus

Development of methodology: D. Bu, C. Lewis, A.F. Gazdar, D.M. Euhus

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): C. Lewis, V. Sarode, A. Lazorowitz, R. Rao, M. Leitch, V. Andrews, A.F. Gazdar, D.M. Euhus

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): D. Bu, C. Lewis, M. Chen, X. Ma, A. Moldrem, A.F. Gazdar, D.M. Euhus

Writing, review, and/or revision of the manuscript: C. Lewis, V. Sarode, R. Rao, V. Andrews, A.F. Gazdar, D.M. Euhus

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): D. Bu, C. Lewis, A. Moldrem, D.M. Euhus

Study supervision: D.M. Euhus

The BRCA1-methylated cell line, UACC3199, was kindly provided by the Arizona Cancer Center, Tucson, Arizona.

This work was financially supported by the Department of Defense Breast Cancer Research Program, contract number W81XWH-07-1-0262, to D. M. Euhus.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.