A New Role for ERα: Silencing via DNA Methylation of Basal, Stem Cell, and EMT Genes

This article is featured in Highlights of This Issue, p. 115

Estrogen receptor α (ERα, ESR1) has proven to be the single most important target in breast cancer. Approximately 70% to 80% of breast cancers are ERα⁺, for which routine testing is used to predict response to antihormonal therapy (1). As demonstrated by genome-wide studies, ERα is a global regulator of gene transcription in breast cancer that orchestrates well-integrated hormonal responses that promote proliferation and survival and inhibit apoptosis (2–5).

As a result of regulating expression of thousands of genes, the presence of ERα drives the luminal classification of breast cancer. There are five intrinsic tumor subtypes, luminal A, luminal B, HER2-enriched, claudin-low, and basal-like, as well as a normal breast-like group. Patients with either luminal B, HER2-enriched, basal-like, or claudin-low tumors experience worse clinical outcome than patients with luminal A tumors (6–8).

ERα has been shown to negatively regulate gene expression, but not much is currently known on how it can achieve this. Epigenetic marks, namely DNA methylation of cytosine at specific CpG sites (5mCpG), are frequently associated with ERα⁺ status in human breast cancers. ERα may play a role in directing DNA methylation to target genes, as specific 5mCpG marks associate with ERα status in human breast cancer and predict risk of tumor recurrence (9–12).

Methylation of cytosine at CpG dinucleotide sites (5mCpG) by DNA methyltransferases (DNMT) in transcriptional regulatory regions mediates stable epigenetic gene silencing. In cancer cells, DNA methylation is highly correlated with repressive chromatin marks, such as trimethylated H3K27 (H3K27me3; ref. 13). H3K27 trimethylation is catalyzed by EZH2, the histone methyltransferase enzymatic subunit of the polycomb repressor complex 2 (PRC2; ref. 14). Together, EZH2 and PRC2 then recruit DNMTs (13, 15). Methylated CpG sites near transcriptional start sites (TSS) can silence gene expression by interacting with effectors, such as methyl-CpG–binding domain proteins, that impede binding of transcription factors, block transcriptional initiation, and recruit histone deacetylases (HDAC) to promote chromatin compaction (16).

The relationship between ERα and DNA methylation patterning in breast cancer has been reported. In a comprehensive bioinformatics study, methylation of CpG sites near ERα-binding regions tended to be lower in ERα⁺ tumors than ERα⁻ tumors. This indicated a passive role for ERα in preventing gene silencing. The methylation status of DNA sequences at ERα-binding sites is tightly coupled with ERα activity (12). Differentially methylated genes have also been identified in antihormone-resistant versus wild-type MCF-7 cells (17, 18), and in ERα RNAi-depleted versus nondepleted MCF-7 cells (19). Consistent with this notion, loss of ERα activity leads to silencing of estrogen-responsive genes, such as PgR (18, 19). Yet, ERα may also play an active role in promoting silencing. A functional link between ERα and DNA hypermethylation has been demonstrated at the CYP1A1 locus, whose gene product converts 17β-estradiol (E₂) into a metabolite that inhibits proliferation; ERα silenced CYP1A1 by recruiting DNMT3B (20).

We sought to identify ERα targets for CpG methylation–mediated silencing by selecting the intersection of (i) genes upregulated (i.e., derepressed) by the demethylating agent decitabine (DAC); (ii) genes upregulated by loss of ERα expression in a series of antihormone-resistant T47D and MCF7 cell lines; and (iii) genes downregulated by E₂ reexposure or increased ERα expression in antihormone-resistant T47D and MCF7 cells. Additional experiments verified the functional dependence on ERα for silencing and DNA methylation of the basal breast cancer subtype markers LCN2 and IFI27 in wild-type and antihormone-resistant T47D-based cell lines. Therefore, we show that ERα targets genes for DNA methylation–mediated silencing that may potentially be predictive biomarkers of breast cancer subtypes.

Sources and culture conditions of cell lines generated in this study are provided in Supplementary Materials and Methods. A schema representing the derivation of antihormone-resistant cells is shown in Supplementary Fig. S1. The fulvestrant (FUL)–resistant cell lines (T47D/FUL, MCF7/FUL) and the estrogen deprivation (ED)–resistant cell lines (T47D/ED1, T47D/ED2) were generated by continuous culture (8 weeks to >1 year) of wild-type T47D and MCF-7 cells in estrogenized media (RPMI1640 plus 10% whole FBS) supplemented with 100 nmol/L fulvestrant or in estrogen-free media (phenol red-free RPMI1640 plus 10% dextran-coated charcoal-stripped FBS), as appropriate. Antihormone-resistant cells were maintained as polyclonal populations. All cell lines were authenticated by gene expression microarrays, morphology, and by verifying ERα, PgR, HER2, LCN2, and IFI27 levels and cell line growth responses to estrogen, ED, and fulvestrant.

The lentiviral cell lines, T47D/ED1/VC, T47D/ED1/VC+E2, T47D/ED1/ERα, and T47D/ED1/ERα+E2 were generated by infecting ERα⁻ T47D/ED1 cells with an ERα-expressing lentivirus or an empty vector control (VC) lentivirus, as appropriate. Infected cells were maintained in estrogen-free or in 1 nmol/L E₂–supplemented medium for 12 weeks. After initial recovery from infection and again 4 weeks later, infected cells were sorted for the lentiviral ZsGreen fluorescent marker using a Becton Dickinson FACS-VantageSE/DiVa cell sorter. To produce the lentiviral vectors, ERα's coding region was excised from pHEGO using EcoRI and inserted into the EcoRI site of the lentiviral vector pLVX-EF1α-IRES-ZsGreen1 (Clontech Laboratories).

RNA was purified using Qiagen's RNeasy Plus Kits. RNA samples were required to exhibit an RNA integrity number of 9.8 to 10.0 on an Agilent 2100 Bioanalyzer.

qRT-PCR assays were carried out as described previously (21) but using AMV First-Strand cDNA Kit, predesigned TaqMan assays, TaqMan Universal PCR Master Mix, and a 7900HT Fast Real-Time PCR system (Thermo Fisher Scientific). Data were analyzed by comparison with a serial dilution series of cDNA. All qPCR data represent the mean and SDs of three independent biological replicates and two technical replicates per biological replicate.

Genome-wide RNA profiling was carried out by the Genomics Facility at Fox Chase Cancer Center (Philadelphia, PA) using Agilent's Human Gene Expression 4 × 44 K v2 oligonucleotide microarrays. RNA labeling (one-color cyanine 3-CTP), hybridization to the arrays, and quality assessment of hybridizations were performed according to the manufacturer's instructions.

Immunoblots were done as described previously (21) but using RIPA buffer and 40 μg protein per lane. Antibodies used are listed in Supplementary Materials and Methods. Blots were visualized using the Odyssey Infrared Imaging System (Li-Cor Biosciences).

Genomic DNA was isolated using the DNeasy Blood and Tissue Kit (Qiagen) and treated with bisulfite [EpiTect Bisulfite Conversion Kit (Qiagen)] to change unmethylated cytosine nucleotides to thymines. Pyrosequencing reactions were carried out at EpigenDx as a service using their predesigned assays. In pyrograms, the ratios of methylated cytosines to thymines (which represent unmethylated cytosines) are internally normalized values. All pyrosequencing data represent the mean and SD of 4 replicates.

Breast cancer data from The Cancer Genome Atlas (TCGA) project were downloaded via the International Cancer Genome Consortium (ICGC) data portal (https://dcc.icgc.org/releases/release_18/Projects/BRCA-US). Methylation data were retrieved for 1,013 patients, 967 of whom also had ERα status available. CpG differential methylation by ERα status was assessed as described in Supplementary Materials and Methods.

The gene expression microarray cohort comprising 2,116 breast cancers was previously constructed in ref. 22 from publicly available datasets and is described in Supplementary Materials and Methods. Breast cancer–intrinsic subtypes were previously classified as described in ref. 22.

To analyze the composite expression level of gene sets in a tumor, gene sets were represented as metagenes and metagene scores, or single number summary values were determined across the expression array breast cancer cohort. These scores represent a linear combination of expression values of each gene in the gene set in individual tumors. Metagene scores were generated by determining the first “principle component” or “eigenvector” of each gene set in each tumor using singular value decomposition (SVD). The eigenvector produced by SVD was rescaled to a rank-based score between 0 and 1, with 0 relating to the lowest composite expression value for a gene set, and 1 relating to the highest. Thus, metagene scores capture the majority of variation in gene expression that is common to the majority of genes in a gene set across a population of samples.

To construct the ERα DNA methylation metagene, Entrez gene identifiers were used to match the Agilent probes from the expression microarrays used in this study to the Affymetrix probe sets used in the combined breast cancer cohort. This resulted in matching 34 of 39 ERα DNA methylation genes (Supplementary Table S1). The ERα status–associated metagene consists of the 100 most differentially expressed genes between ER⁺ and ERα⁻ tumors in the 2,116 breast cancer cohort as determined using the “limma” package (23) for R software (www.r-project.org). The specific genes comprising each metagene are provided in Supplementary Excel File S1.

Microarray data are deposited in the NCBI GEO repository with accession number GSE85536.

Expression array data were log₂ transformed for all comparisons. Differentially expressed genes were identified by serial pairwise comparisons using SAM (24) at an FDR <5% and a 2-fold cutoff, except a 1.5-fold cutoff was used when comparing T47D/ED2/E2 versus T47D/ED2 cells because the ERα levels in these cells were <5% that of wild-type T47D cells (Fig. 1A). Gene enrichment in Supplementary Excel Files S4, S5, and S10 was assessed by one-way Fisher exact tests using the R software application. Associations between ERα DNA methylation metagene scores and distant metastasis–free patient survival (DMFS) were evaluated by Kaplan–Meier analysis and log-rank tests, and by univariate and multivariable Cox proportional hazards regression models as described previously (22). The additional covariates used were age at diagnosis, intrinsic subtype, ERα status, tumor size, and tumor grade. Statistical tests used in Figs. 5 and 6 are specified in the figure legends and were carried out using Prism v4.03 (GraphPad Software). Where specified, one-way tests were employed because gene expression and CpG methylation was assumed a priori to be inversely related.

To identify ERα targets for DNA methylation–mediated silencing, we sought to find the intersection of genes that fulfilled three conditions: (i) those genes derepressed by loss of ERα expression; (ii) those genes resilenced by increased ERα activity or expression; and (iii) those genes derepressed by loss of DNA methylation.

To begin, breast cancer cell line models were developed that exhibited loss of ERα to enable subsequent identification of genes, which inversely correlated with ERα expression/activity. We elected not to use RNAi-based methods, as acute depletion of ERα in estrogen-dependent cells leads to widespread cell death (unpublished observation). Therefore, starting with wild-type ERα⁺ T47D and MCF-7 luminal breast cancer cells, a panel of ERα-low/negative T47D and MCF7 breast cancer cells was derived by long-term selection of cells in 100 nmol/L fulvestrant or in estrogen-free media for 8 weeks to greater than 1 year (schema in Supplementary Fig. S1). Thus fulvestrant-resistant (T47D/FUL, MCF7/FUL) and ED-resistant (T47D/ED1, T47D/ED2) cell lines were derived. ERα mRNA levels were measured by qRT-PCR (Fig. 1A). T47D/ED1 cells lost 99.9%, T47D/FUL and T47D/ED2 cells lost ≥95%, and MCF7/FUL cells (at week 8 of derivation) lost 90% of ERα mRNA compared with respective wild-type parental cells. Immunoblotting also demonstrated similar ERα protein losses (Fig. 4).

To determine global changes in gene expression, which correlated with loss of ERα expression, transcriptional profiling was performed using Agilent 4 × 44 K v2 oligonucleotide microarrays. ERα-low/negative cell lines, T47D/FUL, ED1, ED2, and MCF7/FUL (week 8), were compared against their respective wild-type parental T47D or MCF-7 cell line (four separate pairwise comparisons). Supplementary Excel Files S2 and S3 list the 324 and 153 significantly up- and downregulated genes, respectively. These differentially regulated genes were examined for enrichment of functional gene groups consistent with acquired antihormone resistance using one-way Fisher exact tests (tables and P values in Supplementary Excel Files S4 and S5; Supplementary Fig. S2). As expected, genes with ERα-binding sites were overrepresented. ERE-regulated genes were taken from previously published datasets (see Supplementary Fig. S3 and its legend for a list of genes and references). Importantly, basal markers were very significantly enriched among the upregulated genes, whereas luminal markers were very significantly enriched among the downregulated genes. Also observed was enrichment of cancer stem cell (CSC), epithelial–mesenchymal transition (EMT), and tumor suppressor genes (TSG; see the legend of Supplementary Fig. S3 for references used). Examples of key genes in these functional groups are shown in Supplementary Fig. S3. This indicates that the antihormone-resistant ERα-low/negative T47D and MCF-7 cells transitioned to a differentiation state similar to the basal-like and claudin-low breast cancer subtype. Such a change in differentiation has previously been observed in T47D tumors in vivo following antiestrogen treatment or estrogen withdrawal and termed “luminobasal” (25).

To further refine the list of ERα inversely correlated genes, T47D/ED2 cells were reexposed to E₂ for 38 weeks, resulting in T47D/ED2/E2 cells. Interestingly, ERα RNA (Fig. 1) and protein levels (Fig. 4) never rebounded, indicating permanent ERα silencing as observed elsewhere (26). In fact, ERα RNA levels actually decreased ∼50% more; this likely reflected a known E₂–ERα negative feedback regulatory loop indicative of ERα transcriptional activity (27). MCF/FUL cells were also further selected. These cells at week 8 of derivation showed 90% loss of ERα, but after 13 weeks of additional exposure to fulvestrant (total 21 weeks), ERα levels rebounded to wild-type cell levels (Figs. 1A and 4E). Transcriptional profiling showed increased expression of well-known E₂-stimulated genes in T47D/ED2/E2 versus T47D/ED2 cells (e.g., PGR, CA12, ERBB4) and in MCF7/FUL week 21 versus week 8 cells (e.g., CXCL12, GREB1, ERBB4), as well as decreased expression of E₂-repressed genes (e.g., OASL, C3; both cell lines). Furthermore, the expression pattern of many (but not all) basal and luminal, CSC, EMT, and TSGs reversed upon E₂ reexposure in T47D/ED2/E2 or increased ERα expression in MCF7/FUL (week 21) cells compared with respective parental cells (Supplementary Fig. S3).

Taking into account all cell line transcriptional profiles, 161 genes were identified that consistently inversely related with ERα expression/activity, whereas only 9 genes were directly related (Supplementary Excel Files S6 and S7, respectively).

As ERα inversely related genes whose expression was regulated by DNA methylation were sought, genes upregulated by the DNA demethylating agent DAC were identified. Wild-type T47D cells were treated with 1 μmol/L DAC or control (CON)-treated for 96 hours and then transcriptionally profiled. This resulted in the identification of 1,049 genes (Supplementary Excel File S8).

Subsequently, the intersection of ERα inversely related genes and DAC-induced genes was determined. This intersection represented the set of genes that fulfilled the following criteria: (i) genes induced by DAC versus CON-treated wild-type T47D cells; (ii) genes upregulated in each of the ERα-low/negative cell lines, that is, T47D/FUL, T47D/ED1, and T47D/ED2 cells, versus wild-type T47D cells; (iii) genes downregulated by E₂ in T47D/ED2/E2 versus T47D/ED2 cells; (iv) genes upregulated in ERα-low MCF7/FUL week 8 versus wild-type MCF-7 cells; and (v) genes downregulated in ERα⁺ MCF-7/FUL week 21 versus ERα-low MCF-7/FUL week 8 cells. These selection criteria pinpointed 39 high-value candidates for ERα-mediated silencing via DNA methylation (Fig. 1B; Supplementary Excel File S9).

Initially, these 39 candidate genes were evaluated for methylation in human breast cancer (Supplementary Table S1). Using TCGA-processed breast cancer methylation data, a set of 1,996 CpGs associated with these genes was identified. These CpG sites were assessed for differential methylation between ERα⁺ and ERα⁻ breast cancers using one-sided Wilcoxon rank sum tests adjusted for FDR. Using a permutation analysis to determine whether similar results could be achieved using 1,000 sets of 39 random genes, it was concluded that the candidate ERα DNA methylation targets tended to display higher methylation levels in ERα⁺ compared with ERα⁻ tumors than would be expected for an identically sized set of randomly selected genes (permutation-based P = 0.011).

Next, the 39 candidate ERα DNA methylation targets were analyzed for enrichment of the same gene groups as the ERα inversely related genes (tables and P values in Supplementary Excel File S10; Supplementary Fig. S2). Similar to the earlier results, the candidate methylation targets were enriched for genes with ERα-binding sites, basal markers, CSC-upregulated genes, EMT-upregulated genes, and TSGs. The candidate methylation targets were also enriched for EMT-downregulated genes, but there were almost twice as many EMT-upregulated genes than downregulated (13 vs. 7, respectively).

The candidate ERα DNA methylation target gene set was analyzed relative to other important tumor-related gene sets and clinical variables in a cohort of 2,116 breast cancers. Gene sets were represented as a composite entity termed an “expression metagene” and a single value summary of the gene set's expression level in an individual tumor as a “metagene score.” To enable evaluation of the distribution of gene set expression levels, metagene scores were used to divide the breast cancers in the cohort into “tertiles” (lowest 33%, middle 33%, highest 33%).

ERα DNA methylation metagene scores were plotted versus ERα status, an ERα status–associated metagene, breast cancer intrinsic subtypes, luminobasal signature metagenes, EMT metagenes and CD44⁺/CD24^−/low CSC metagenes (Fig. 2). The ERα status–associated metagene encapsulated the 100 most differentially expressed probe sets between ERα-positive and -negative tumors in the 2,116 breast cancer cohort. Congruent with the original selection criteria, the ERα DNA methylation metagene showed a clear negative association with ERα status and 100 other ERα status–associated genes, with ERα⁺ tumors tending to have lower scores (and thus indicating lower levels of gene expression; Fig. 2A). With regard to intrinsic subtype, luminal A and B subtypes displayed the lowest ERα DNA methylation metagene scores, whereas the basal-like subtype exhibited the highest scores (Fig. 2B). This was consistent with enrichment of basal-up/luminal-down genes as previously noted. Furthermore, the ERα DNA methylation metagene clearly directly related to the luminobasal signature metagenes (Fig. 2C), suggesting the ERα DNA methylation targets program this type of change in differentiation. Again as expected from the enrichment analysis, ERα DNA methylation metagene scores were associated with CD44⁺/CD24^−/low metagenes (Fig. 2D) and selectively with the EMT-upregulated metagene (Fig. 2E). This helps explain why the ERα DNA methylation metagene scores also associated with the claudin-low breast cancer subtype (Fig. 2B).

The ERα DNA methylation metagene was next evaluated for predicting DMFS in the breast cancer cohort by Kaplan–Meier survival curves (Fig. 2F) and Cox proportional hazards regression models (Supplementary Table S4). In each analysis, patients were separated into two groups according to metagene scores split at the 10th, 50th, or 90th percentiles, and then the proportion of patients exhibiting DMFS in each group was plotted against time. Patients with metagene scores in the bottom 10th or top 90th percentile experienced significantly decreased DMFS. Likewise, univariate Cox proportional hazards regression analysis demonstrated that metagene scores split at the 10th and 90th percentiles associated with DMFS (P = 0.00003 and 0.035, respectively). ERα DNA methylation metagene scores split at the 10th percentile remained significantly associated with DMFS in a multivariable Cox proportional hazards regression model (P = 0.026), but not when split at the 90th percentile. These results suggested that some genes in the ERα DNA methylation metagene when expressed at low levels promoted poor DMFS, whereas others did so when expressed at high levels.

To determine which of the genes of the candidate ERα DNA methylation metagene when expressed at low or high levels may promote poor DMFS, Kaplan–Meier survival curves and univariate Cox proportional hazards regression analysis was conducted for each gene. On the basis of these analyses, the candidate ERα methylation targets were separated into low- and high-expression metagenes (defined in Supplementary Excel File S11). The low- and high-expression metagenes poorly correlated with each other, indicating they indeed likely represented different biological processes (Supplementary Table S5). Patients were then divided according to their tumor's low- and high-expression metagene scores split at the 50th percentile and evaluated for DMFS as before (Fig. 3A). The Kaplan–Meier plots showed clear separations of survival curves in which patients in the low expression metagene's bottom 50% group and patients in the high expression metagene's top 50% group displayed poor DMFS.

Low- (Fig. 3B) and high-expression ERα DNA methylation metagenes (Fig. 3C) were next assessed for associations with various tumor-related metagenes in the breast cancer cohort. These results indicated that the low-expression metagene associated with tumor suppressor and focal adhesion gene expression in breast cancers. Accordingly, low levels of these types of genes would be predicted to promote metastasis. Conversely, the high-expression metagene associated with high-grade tumors, as well as proliferation and proinflammatory Th1 immune response gene expression in breast cancers; this would also promote poor DMFS. Supplementary Excel File S11 contains references that help provide a rationale for the segregation of genes into either the low- or high-expression metagenes.

IFI27 and LCN2 were the top two genes inversely related to ERα expression/activity in the T47D-based cell lines (Fig. 1B). Both LCN2 (5) and IFI27 (2) contain ERα-binding sites. Also, both are basal markers (7, 28, 29) and promote EMT (30, 31). Hence, LCN2 and IFI27 were selected for validation of ERα-dependent changes in expression and 5mCpG levels.

LCN2 mRNA and protein levels dramatically increased in a time-dependent manner after precipitous drops in ERα levels across all antihormone-resistant models (Fig. 4A–C, and E). Furthermore, LCN2 expression decreased in a time-dependent manner after extended E₂ reexposure in T47D/ED2/E2 cells (Fig. 4D) and once ERα expression rebounded in MCF/FUL cells (Fig. 4E). LCN2 has previously been reported to upregulate the key EMT transcription factor slug (SNAI2; ref. 30); therefore, slug expression was examined. Across all antihormone-resistant models, changes in slug expression followed similar changes in LCN2, although in MCF7/FUL cells, slug induction was delayed until after LCN2 was upregulated from 8 to 12 weeks and silenced again due to ERα reexpression (Fig. 4B–E). Together, these results are consistent with LCN2 regulating slug expression and demonstrate that LCN2 may have promoted EMT via slug.

Like LCN2, IFI27 mRNA expression was strikingly upregulated 270- to 1,900-fold across the ERα-low/negative T47D-based antihormone-resistant compared with wild-type cell lines and was dramatically repressed again upon E₂ reexposure in T47D/ED2/E2 versus T47D/ED2 cells (Fig. 4F).

Levels of selected 5mCpG sites near the TSSs of LCN2 and IFI27 were quantitated across the T47D-based models by pyrosequencing of bisulfite-treated gDNA. This analysis found LCN2 and IFI27 CpG sites to be significantly hypomethylated in ERα-low/negative T47D/FUL, T47D/ED1, and T47D/ED2 cells versus wild-type ERα⁺ T47D cells (Fig. 5A), and significantly hypermethylated upon E₂ reexposure in T47D/ED2/E2 cells compared with parental T47D/ED2 cells (Fig. 5B). Therefore, LCN2 and IFI27 CpG methylation levels directly associated with ERα expression/activity.

Next, a causal relationship between ERα expression and CpG methylation of LCN2 and IFI27 was tested. ERα⁻ T47D/ED1 cells were infected with an ERα-expressing lentivirus or an empty VC lentivirus generating T47D/ED1/ERα and T47D/ED1/VC cells, respectively. These infected cells were maintained with and without E₂ for 12 weeks and subjected to two rounds of cell sorting for the lentiviral ZsGreen fluorescent marker. Characterization of these lentiviral cells lines demonstrated functional ERα signaling as ERα and PgR mRNAs were expressed at high levels in ERα-infected cells. Also, LCN2 and IFI27 mRNA levels were downregulated in an ERα-dependent manner (plus E₂ for IFI27; Fig. 5C). Next, LCN2 and IFI27 CpG methylation levels were quantitated by pyrosequencing and found to be significantly increased in lentiviral ERα plus E₂ (LCN2) or just ERα (IFI27) compared with VC cells. (Fig. 5D). Therefore, increased CpG methylation of LCN2 and IFI27 was dependent on ERα plus E₂ stimulation. In the case of IFI27, repression of its expression did not occur until its CpG methylation levels were maximally increased by the presence of E_2, indicating that perhaps a methylation threshold was needed to cause its repression.

LCN2 and IFI27 were examined in a panel of 11 breast cancer cell lines. Initial characterization showed 4 cell lines were ERα⁺ and 4 were HER2⁺, where only BT-474 cells were positive for both prognosticators (Fig. 6A). Then, LCN2 protein (Fig. 6B) and IFI27 mRNA levels (Fig. 6C) were measured and found to be significantly lower in ERα-positive compared with ERα⁻ cells. Next, methylation levels of CpG sites near the TSSs of LCN2 and IFI27 were quantitated by pyrosequencing of bisulfite-treated gDNA. Correlations between expression and 5mCpG levels were determined by Spearman ρ (Fig. 6D). Methylation of all 5 of LCN2's tested CpG sites and 4 of IFI27′s 9 tested sites (CpG sites +404, +438, +508, and +550) significantly correlated with each respective gene's expression levels. Finally, an association between CpG methylation levels and ERα status was evaluated. Although all 5 of LCN2′s CpGs were evaluated, only those 4 CpGs of IFI27 that significantly correlated with expression were considered. This analysis showed for both LCN2 and IFI27, that 5mCpG levels were significantly higher in ERα⁺ than ERα⁻ cells (Fig. 6E).

LCN2 and IFI27 RNA expression levels were examined in the breast cancer cohort with respect to DMFS by Cox proportional hazards regression analysis (Supplementary Table S6). In univariate models, LCN and IFI27 both significantly associated with DMFS (P = 0.040 and 0.0023, respectively), but this did not hold in multivariable models.

We hypothesized that ERα may regulate gene expression in part via DNA methylation as methylation of specific CpG sites associates with ERα⁺ status in human breast cancer. This hypothesis was tested by identifying genes normally silenced in ERα⁺ breast cancer cell lines but which were derepressed upon exposure to the demethylating agent DAC, derepressed upon long-term loss of ERα expression, and resuppressed by gain of ERα activity/expression. On the basis of these criteria, 39 candidate ERα DNA methylation targets were found. These 39 targets were used to construct an ERα DNA methylation metagene that inversely associated with ERα status in human breast cancers and directly associated with expression signatures of basal-like and claudin-low breast cancer subtypes (25). Congruent with these associations, the candidate ERα DNA methylation targets were enriched for basal markers, CSC, and EMT genes.

LCN2 and IFI27 were the top two ERα inversely related genes identified and were selected for validation. Both LCN2 (5) and IFI27 contain ERα-binding sites (2), are basal markers (7, 28, 29), and involved in EMT (30, 31). First, LCN2 and IFI27 were originally silenced in wild-type T47D and MCF-7 cells, but their expression dramatically increased upon loss of ERα, while their 5mCpG levels significantly decreased in all antihormone-resistant T47D cell lines. Second, LCN2 and IFI27 were resilenced upon E₂ reexposure in T47D/ED2 cells while their 5mCpG levels increased. Third, lentiviral ERα plus E₂ in T47D/ED1 cells also repressed LCN2 and IFI27 expression while increasing their 5mCpG levels. Fourth, LCN2′s and IFI27′s 5mCpG levels positively associated with ERα status but inversely correlated with expression in a panel of 11 breast cancer cell lines. Together, these results provide correlative and functional evidence that ERα directed DNA methylation–mediated silencing of LCN2 and IFI27.

As ERα plays such a pivotal role in a more favorable outcome in breast cancer, genes targeted by ERα for DNA methylation–mediated silencing likely play important roles in disease progression. In addition to the CSC and EMT genes, Kaplan–Meier survival curve analyses indicated that the candidate ERα DNA methylation targets consisted of two classes of genes that predicted poor DMFS, one when expressed at low levels and a second when expressed at high levels. The low expression class associated with tumor suppressor and focal adhesion gene expression in breast cancer. Conversely, the high-expression class associated with proliferation and inflammatory response gene expression in breast cancer. In addition, the two validated targets for methylation, LCN2 and IFI27, predicted DMFS in univariate Cox proportional hazards models. Thus, genes methylated and silenced in an ERα-dependent manner may be good targets for therapeutic intervention in ERα⁻ breast cancer where they are expressed.

How might ERα direct DNA methylation to specific genes? We propose it may begin with transcriptional repression. E₂ actually represses transcription of more genes than it stimulates (3, 4). Studies on E₂-dependent transcriptional repression have demonstrated that ERα recruits coregulators (corepressors; refs. 32–36) and coactivators that act as corepressors (36, 37). The coregulators serve as scaffolds to interact with HDACs and a host of additional cofactors (32, 33, 38, 39), such as EZH2 (34, 40), that together remove activating histone marks, add repressive marks, and restructure chromatin structure (38, 39, 41–43). We hypothesize that not only does ERα direct epigenetic silencing via histone modification, but also via cytosine methylation at CpG sites.

Genome-wide kinetics of DAC-induced DNA demethylation and subsequent remethylation after drug withdrawal in breast cancer cells showed that CpGs differ in both their susceptibility to demethylation and propensity for remethylation after drug removal (44).

This is also plausible as EZH2 recruits DNMTs directly and indirectly through PRC2 (13, 15). Other protein–protein interactions exist as well that could support formation of a multicomponent complex containing ERα and DNMTs, such as those between ERα and EZH2 (45, 46), between HDACs and DNMTs (47–49), and between corepressors and DNMTs (40, 49). Evidence for such a complex exists at least at the CYP1A1 promoter, where it was demonstrated that ERα and DNMT3B interacted (20). Thus, it is possible that ERα could silence targeted genes via DNA methylation by directly and indirectly recruiting corepressors, HDACs, EZH2 in PRC2, and DNMTs (model shown in Fig. 7).

Taken together, our data indicate that ERα can silence genes via DNA methylation, such as LCN2 and IFI27. Moreover, ERα may direct DNA methylation–mediated silencing of a subpopulation of basal markers, CSC, and EMT genes that may potentially enforce luminal differentiation of breast cancer cells.

No potential conflicts of interest were disclosed.

Conception and design: E.A. Ariazi, J. Boyd

Development of methodology: E.A. Ariazi, J.C. Taylor

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): E.A. Ariazi, J.C. Taylor, E. Nicolas, J. Boyd

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): E.A. Ariazi, M.A. Black, M.J. Slifker

Writing, review, and/or revision of the manuscript: E.A. Ariazi, M.A. Black, M.J. Slifker, D.J. Azzam, J. Boyd

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): E.A. Ariazi, J.C. Taylor, M.J. Slifker

Study supervision: E.A. Ariazi, J. Boyd

The authors thank Dennis DeSimone and Trung Nguyen, clinical fellows in the laboratory for technical support. The authors also thank the Expression Microarray facility, the Genotyping and Real-Time PCR facility, and the Flow Cytometry facility at Fox Chase Cancer Center for technical support.

This work was supported by the Commonwealth Universal Research Enhancement (CURE) Program Award from the Pennsylvania Department of Health (to J. Boyd) and NIHP30 CA006927 (Fox Chase Cancer Center Core Grant).

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.