Activation of the Glucocorticoid Receptor Is Associated with Poor Prognosis in Estrogen Receptor-Negative Breast Cancer

Breast cancer, the most frequent cancer in women, is generally classified into estrogen receptor–positive (ER⁺) and estrogen receptor–negative (ER⁻) subtypes because of the strong prognostic and predictive importance of ER-alpha (1, 2). ER-alpha expression (henceforth, ER⁺) is associated with a significantly lower 5-year recurrence rate in part because of the effectiveness of ER-targeted treatments (3) and in part because of a generally more indolent and well-differentiated phenotype when compared with ER⁻ tumors (2, 4). As a consequence, identifying ER-positivity (usually through immunohistochemistry or IHC) has become routine clinical practice. Unfortunately, the progress made in treating ER⁻ breast cancer has been less successful. First, ER⁻ tumors are often intrinsically more aggressive and of higher grade than ER⁺ tumors; second, underlying signaling pathways driving ER⁻ breast cancer have been difficult to identify in part because of ER⁻ tumor heterogeneity. Compared with ER⁺ tumors, relatively few prevention and treatment strategies are available for ER⁻ breast cancer patients (5, 6). Until recently, the only clinically effective treatment for ER⁻ breast cancer was cytotoxic chemotherapy. More recently, antiangiogenic and PARP inhibitors have been explored as treatments for ER⁻ breast cancer (7), although their ultimate effectiveness in the ER⁻ breast cancer subgroup is not yet clear. Despite the fact that many “triple-negative” [i.e., ER⁻, progesterone receptor–negative (PR⁻), and HER2⁻] breast cancers are highly sensitive to chemotherapy, a significant proportion show resistance to chemotherapy and progress rapidly (8).

Recently, the stress hormone (glucocorticoid) receptor (GR) was shown to play an important role in breast epithelial and breast cancer cell biology, particularly in ER⁻ premalignant and breast cancer cell lines (9–11). Using ER⁻/GR⁺ breast epithelial and cancer cell lines (e.g., MCF10A-Myc and MDA-MB-231) and global gene expression studies, our laboratory identified several GR-regulated target genes whose protein products mediate cell survival, including serine/threonine-protein kinase 1 (SGK1; ref. 12) and the dual specificity protein phosphatase 1 (DUSP1/MKP1; ref. 13). In addition, using a mouse xenograft model of MDA-MB-231 breast cancer cell tumor growth, we found that tumor cell apoptosis induced by paclitaxel is inhibited by pretreatment with systemic dexamethasone (dex, a synthetic glucocorticoid; ref. 14). We also discovered that increased endogenous glucocorticoid exposure is associated with greater ER⁻ mammary tumor growth in the C3(1) SV40 Tag transgenic mouse model of ER⁻/PR⁻/GR⁺ invasive breast cancer (15). These observations suggest that GR signaling may be an important pathway influencing ER⁻ breast cancer biology.

To further understand the role of GR signaling in ER⁻ breast cancer, we carried out genome-wide GR chromatin immunoprecipitation-sequencing (ChIP-seq) and time-course gene expression analysis, using the ER⁻/GR⁺ premalignant breast epithelial cell line MCF10A-Myc. We found that GR target genes were enriched in many cancer-related pathways, including epithelial-to-mesenchymal transition (EMT), cell adhesion, and cell survival. In addition, to examine whether GR-mediated gene expression associates with worse outcome in ER⁻ breast cancer patients, we compiled and analyzed a meta-dataset of early-stage breast cancer patients with associated tumor gene expression and long-term clinical outcome data. In this retrospective analysis of primary tumor gene expression data in 1,378 patients, we found that high GR gene (NR3C1) expression was associated with a significantly shorter relapse-free survival (RFS) in both adjuvant chemotherapy–treated and untreated early-stage ER⁻ patients. Unexpectedly, we also found that in ER⁺ breast cancer, higher GR expression was associated with longer RFS. Thus, both the cell line and primary human tumor data analyses suggest that the GR directly mediates gene expression pathways associated with aggressive behavior and early relapse in the ER⁻ breast cancer subtype.

Early passage MCF10A-Myc cells were derived from MCF10A cells (ATCC), transduced with the viral oncogene c-Myc as previously described (9), and checked for Myc expression. Cells were cultured in a 1:1 mixture of Dulbecco's modified Eagle's medium and Ham's F12 (BioWhittaker) supplemented with hydrocortisone (0.5 μg/mL), human recombinant epidermal growth factor (10 ng/mL), insulin (5 ng/mL), 5% FBS, 100 units/mL penicillin, and 100 μg/mL streptomycin. Identical conditions were used for gene expression profiling after dex treatment (16).

MCF10A-Myc cells (1 × 10⁷) were serum starved for 3 days and then treated with either dex or ethanol for 1 hour as previously described (9). Cellular DNA and protein were then cross-linked for 20 minutes with formaldehyde (1% final concentration) followed by addition of glycine to a final concentration of 125 mmol/L for 5 minutes. We also took 1% of each lysate to evaluate as input protein and DNA. ChIP experiments were conducted according to a standard protocol (Upstate Biotechnology, Millipore). ChIP-grade rabbit polyclonal anti-GR (sc-1003x; Santa Cruz) antibodies were used for immunoprecipitation. We also used normal rabbit IgG (sc-2027; Santa Cruz) as a negative control.

An anti-GR Western blot analysis was carried out on a small percentage of the input lysate and the anti-GR immunoprecipitate to assess enrichment of GR by the immunoprecipitation process (Supplementary Fig. S1). Cell lysates or immunoprecipitated complexes were resuspended in 2× Laemmli lysis buffer. The proteins were resolved by SDS-PAGE, transferred to nitrocellulose membrane, blocked with 5% skim milk in TBS-T (0.1% Tween-20 in TBS), and probed with a 1:1,000 dilution anti-GR antibody (sc-1003; Santa Cruz). Goat anti-rabbit horseradish peroxidase (HRP; Santa Cruz; 1:5,000) was used as the secondary antibody. Proteins were visualized by enhanced chemiluminescence (Amersham Pharmacia Biotech). After confirming efficient GR immunoprecipitation (Supplementary Fig. S1A and B), we reversed the cross-links of the immunoprecipitated protein–chromatin complexes according to the manufacturer's instructions. ChIP DNA was then purified by Qiagen's PCR Purification KIT.

To ensure the quality of the anti-GR antibody ChIP'd DNA, we carried out quantitative real-time PCR (qRT-PCR) as previously described (13) on the promoter regions of SGK1 and TSC22D3/GILZ genes containing established GR binding regions (GBR; refs. 17, 18). The primers for SGK1′s promoter regions were 5′-CCCCTCCCTTCGCTTGTT-3′ (sense) and 5′-GGAAGAAGTACAATCTGCATTTCACT-3′ (antisense; ref. 19); the primers for TSC22D3/GILZ's promoter regions were 5′-GATACCAGTTAAGCTCCTGA-3′ (sense) and 5′-AGGTGGGAGACAATAATGAT-3′ (antisense; ref. 20). Dex- and ethanol-treated ChIP DNA and input samples were analyzed in triplicate.

After the quality control experiments described earlier proved to be successful, the GR ChIP DNA was deep-sequenced by using Illumina's Solexa sequencer [National Center for Genome Resources (NCGR)]. We used the Maq program (version 0.71; ref. 21) to align the resulting sequenced 36-bp tags to the Human Genome (NCBI/b36). Tags that were mapped to more than one location in the human genome were removed. The GBRs were identified using the Model-based analysis of ChIP-Seq (MAC) program (version 1.3.7.1; ref. 22). Input DNA tags were used to call GBRs.

A P value cutoff of 10⁻³ or less was initially used to call the GBRs in both dex- and ethanol-treated samples. We then developed more stringent criteria for accurately detecting significant peak signals based on known GR promoter occupancy regions for the 2 GR target genes, SGK1 (19) and TSC22D3/GILZ (20). In the initial P ≤ 10⁻³ GBR dataset, we found the number of tags per binding region (N_tag) in all the SGK1 and TSC22D3/GILZ–associated GBRs in ethanol-treated samples was always 6 or less, whereas N_tag in dex-treated samples was always 5 or more (see Supplementary Table S1 for the summary of SGK1 and TSC22D3/GILZ GBR data). Because ethanol-treated samples should not show significant GR occupancy for these 2 previously described dex-dependent GBRs, an N_tag > 6 was determined as the cutoff for removing non-dex-dependent GBRs associated with ethanol treatment. Similarly, we identified −log₁₀(P) > 5.8 [the maximum value of −log₁₀(P) in ethanol-treated samples] as the minimum cutoff for removing likely nonspecific GBRs associated with ethanol treatment. Therefore, N_tag > 6 and −log₁₀(P) > 5.8 became our lower-bound cutoffs for removing likely false-positive (i.e., ethanol-associated) GBRs in the MCF10A-Myc cell ChIP-seq experiments.

We also validated the dex-dependent GR occupancy of the same established promoter region GBRs in SGK1 and GILZ, using a directed ChIP-PCR assay (Supplementary Fig. S1B). In the same sample subjected to dex-treated ChIP-seq, we found that the N_tag for these 2 GBRs were 22 (SGK1) and 25 (GILZ). For the promoter GBR for SGK1 and GILZ, the −log₁₀(P) were 13.9 and 15.5, respectively. Both values were much higher than the lower-bound cutoffs used to filter out ethanol-associated GBRs (i.e., N_tag > 6 and −log₁₀(P) > 5.8). Because SGK1 and GILZ GBRs are known to have very high GR occupancy, N_tag > 22 and −log₁₀(P) > 13.9 were likely overly stringent cutoffs for genome-wide identification of GBRs. We therefore took the average of the N_tag for the lower-bound cutoff and for the more stringent cutoff as the final cutoff [(6 + 22)/2 = 14], i.e., N_tag ≥ 15]. Likewise we averaged the −log₁₀(P) cutoffs to determine the final cutoff value [(5.8 + 13.9)/2 = 9.85 and rounded to 10], that is, −log₁₀(P) ≥ 10 or P ≤ 10⁻¹⁰. The final set of GBRs was then identified by using both N_tag ≥ 15 and P ≤ 10⁻¹⁰.

All other statistical analyses and data plotting were carried out using the R language. Gene pathways were identified by MetaCore database and software suite version 6.6 build 28323 (GeneGo Inc.) Enrichment analysis of transcription factor binding motifs located in GBRs was carried out by TRANSFAC version 2009.3 (23, 24).

Microarray raw data CEL files of the MCF10A-Myc cell time-course (2, 4, and 24 hours) gene expression profiling following dex or ethanol treatment (16) were downloaded from the Gene Expression Omnibus (GEO) database (GEO ID is GSE4917). We then renormalized and reanalyzed the original data by using the latest version of Robust Multiarray Average (RMA) method (25) in the Bioconductor's affy package (26). Up- and downregulated genes were identified as those with greater than a ±1.5 fold change in the average of dex- or ethanol-treated values from 3 independent experiments at any one time point.

Breast cancer gene expression raw data CEL files derived from early-stage primary tumors with clinical follow-up from 8 well-annotated studies (Table 1) were downloaded from GEO. To minimize potential cross-platform bias, we used data obtained only from Affymetrix's U133A and/or U133+2 platforms (Affymetrix's U133A and U133+2 are essentially the same chip except that U133+2 contains more probesets than U133A). Duplicate patient samples were identified by both GEO ID and gene expression and removed from the dataset. Assuming that tumor microarray analyses within the same “batch” or hospital/lab source were carried out under similar conditions, we first grouped samples by institution. We then used the normalized unscaled standard error (NUSE) analysis in Bioconductor's affyPLM package (27, 28) to remove poor quality arrays within each institutional batch. Next, we carried out an RMA normalization (25) within each study group, using Bioconductor's affy package (26) to adjust the background noise. To remove batch effects between study groups, we carried out a within-array normalization by subtracting the array average expression of 5 well-established housekeeping genes (ACTB, GAPDH, GUSB, RPLP0, and TFRC, altogether 19 probesets) from the expression of each probeset. We confirmed that the variation of the mean expression of these 5 housekeeping genes was very small in individual datasets (all coefficients of variation were less than 0.03), ensuring that the normalized values were comparable across the groups. A similar normalization strategy for qRT-PCR assessment of mRNA has been used on the commercially available Oncotype DX assay (29).

We carried out a receiver operating characteristic (ROC) analysis by using the Bioconductor's genefilter package (30) on tumors with known ER (IHC) status in the meta-dataset (799 ER⁺ and 208 ER⁻) and confirmed that the ESR1 205225_at probeset was most accurate for predicting ER IHC status (ref. 31; Supplementary Fig. S2A). We then calculated the Youden Index [J = MAX(specificity + sensitivity − 1)] of the 205225_at probeset ROC curve (32, 33), and took the corresponding normalized expression of the probe as the cutoff point to determine ER status. This cutoff value was −3.416, meaning that if ESR1 expression was greater than −3.416, the tumor was identified as ESR1⁺; otherwise it was considered to be ESR1⁻.

Because there are no available GR IHC datasets coupled to gene expression data, we were unable to identify a GR gene (NR3C1) expression cutoff by using ROC analysis. There are 4 NR3C1 probesets: 201865_x_at, 201866_s_at, 211671_s_at, and 216321_s_at on the Affymetrix's U133A and U133+2 arrays. We used 216321_s_at to represent the expression of NR3C1 because (a) only 216321_s_at hybridizes with all NR3C1 mRNA isoforms (Supplementary Fig. S2B); (b) only the probe sequences (11 for each probeset) of 216321_s_at and 211671_s_at can be all aligned perfectly (100% score) with NR3C1 [by BLAT (34), human genome NCBI36/hg18]; (c) on the basis of BioGPS (35), 216321_s_at shows the most diverse tissue expression of the 4 probesets. We then used the highest and lowest quartiles (25%) of NR3C1 probeset expression as a cutoff to identify “high” and “low” GR (NR3C1) tumor expressions, respectively.

RFS was estimated by using the Kaplan–Meier method. RFS between groups was compared by using the log-rank test. Hazard ratios were calculated by using the Cox proportional hazards regression model. All survival analyses were carried out by the survival package in R (36).

Genes that were differentially expressed between NR3C1 (GR)-high and NR3C1 (GR)-low tumors were identified by using Significance Analysis of Microarray (SAM; ref. 37) in R's siggene package (38). We also used the MetaCore software suite (GeneGo Inc.) to identify significant gene pathways from differentially expressed genes between GR-high and GR-low tumors. All other statistical analyses were done using the R language.

To examine the genome-wide GBRs in a model of ER⁻ breast cancer, MCF10A-Myc cells were subjected to GR ChIP-seq (39). The genomic locations of the sequence tags obtained from deep-sequencing were mapped by using the UCSC's human genome version NCBI36/hg18. The number of uniquely mapped 36-bp tags was 1.34 × 10⁷ in the dex-treated sample and 0.83 × 10⁷ in the ethanol-treated sample, suggesting significantly more GR binding events following dex treatment. Using our customized threshold values developed from well-validated GBRs (SGK1 and GILZ), we identified 1,533 GBRs following dex treatment and 126 GBRs following ethanol treatment. Eighteen dex-associated GBRs were found within 100 kb of any ethanol-associated GBRs. The remaining 1,515 GBRs seemed to be highly specific for ligand-bound GR (a list of dex-specific GBRs is provided in Supplementary Table S2). Figure 1 shows that more than half of these dex-specific 1,515 GBRs were located in intergenic regions. Even within genes, most of the GBRs were located in introns, showing the typical distribution of genome-wide binding sites reported previously for ligand-bound nuclear receptors such as ER-alpha (40). GR response elements were found in 66% (1,000 of 1,515) of the dex-dependent GBRs identified, making the canonical GR binding motifs the third most commonly identified transcription factor. The 10 most common transcription factor motifs identified by TRANSFAC analysis are shown in Supplementary Table S3.

GR target gene identification was based on the presence of at least one GBR within 100 kb of a gene's transcription start site (TSS). GR target genes were further classified into 2 categories, based on the physical distance of the GBR from the TSS: (a) 0 to 10 kb (a proximal and traditional “promoter” location) and (b) 10 to 100 kb (a distal and traditional “enhancer” location). Using these criteria, 286 target genes had proximal GBRs and 2,044 genes had distal GBRs; 133 genes showed GR occupancy in both regions. Therefore, GR binding occurred more frequently in the distal/enhancer regions of target genes, although the density of GR occupancy was higher in the more proximal regions [28.6 (proximal) vs. 22.7 (distal) hits/kb].

Because GR chromatin occupancy does not guarantee a transcriptional event, we refined our list of likely direct GR target genes, using time-course gene expression data from MCF10A-Myc cells (GEO accession ID is GSE4917; ref. 16). This represents the steady-state mRNA expression following treatment with dex (10⁻⁶ mol/L) versus ethanol for 2, 4, and 24 hours. In total, 680 genes were significantly regulated (an average of ≥ 1.5 or ≤ −1.5 fold change compared with the ethanol vehicle at any one or more time points). Among them, 517 genes were upregulated and 163 genes were downregulated. There was no overlap between the gene sets of upregulated and downregulated genes.

By examining the MCF10A-Myc cell line ChIP-seq target genes and the gene profiling list, we identified 187 overlapping genes, among which 51 had proximal GBRs, 160 had distal GBRs, and 24 had both. Putative direct target genes identified by overlapping the ChIP-seq–identified genes with the gene profiling lists (either up or downregulated) were compared with the total number of differentially expressed genes at each time point (Fig. 2). This revealed a much larger proportion of directly upregulated genes composing the more immediate 2- and 4-hour time points compared with the later 24-hour time point (Fig. 2A). This finding is consistent with the direct GR-mediated transcriptional activation of these earlier genes and was true regardless of the proximal versus distal location of the GBR. Interestingly, a much lower proportion of repressed genes (Fig. 2B) was identified as direct GR target genes. This suggests that GR more commonly mediates transcriptional repression through indirect mechanisms (e.g., due to sequestration of activating transcription factors such as NF-κB), rather than through direct GR occupancy of “repressive” DNA elements. In the downregulated genes that were identified by gene expression studies and also showed GR occupancy by ChIP-seq, the GR bound almost exclusively to regions greater than 10 kb distal to the TSS, consistent with previous observations in lung cancer cells (41). Finally, unlike upregulated genes which mainly showed significantly increased expression 2 and 4 hours following dex treatment, the proportion of genes with downregulated mRNA steady-state levels following GR activation was approximately equal at all time points, suggesting that the kinetics of GR-mediated gene repression may be significantly influenced by variable mRNA degradation rates following repression of gene transcription.

We next carried out a gene enrichment pathway analysis to identify the functional pathways associated with the 187 GR “direct” target genes that were identified in the MCF10-A Myc premalignant breast cancer cell line following GR activation (Table 2). The top 10 most significant pathways included GnRH activation, EMT, immune regulation, and proliferation, all of which are critical pathways associated with clinical outcome in breast cancer. Taken together, the MCF10A-Myc cell line data suggested that GR activation directly regulates genes associated with aggressive breast cancer biology.

Given the pathways revealed by GR transcriptome analysis above, and because GR activation in ER⁻ breast cancer cell lines has been shown to promote cell survival, chemotherapy resistance, and increased tumor growth in a preclinical xenograft model of human breast cancer (14), we next asked whether high versus low GR expression in primary breast cancers is associated with worse outcome in early-stage ER⁻ breast cancer patients. To test this hypothesis, we compiled a meta-dataset of the 8 publically available Affymetrix gene expression studies of early-stage primary breast tumors in patients with long-term clinical follow-up (refs. 42–49; n = 1,378, summarized in Table 1).

Because the longer time to relapse in ER⁺ versus ER⁻ breast cancers suggests clinically different disease entities, we wished to analyze these tumor subsets independently. However, ER IHC information was not available for all tumors; therefore, we used ER-alpha gene (ESR1) expression to determine ER status. Using the expression of ESR1 probeset 205225_at as determined by ROC analysis (area under the curve = 0.939), we identified 1,024 ESR1⁺ and 354 ESR1⁻ tumors. For GR expression, we were not able to identify a GR gene (NR3C1) expression cutoff using ROC analysis because there are no publicly available GR IHC datasets coupled with NR3C1 gene expression data. Published GR IHC studies, however, have shown a range of 15% to 40% GR-positivity in invasive breast cancers; this variation is likely dependent on tumor subtype, patient age, ancestry, sample conditions, and antibody (10, 50, 51). As described in the Materials and Methods, because most studies show that at least 25% of breast cancers express significant GR by IHC, we used the highest and lowest quartiles (top and bottom 25%) of NR3C1 probeset 216321_s_at expression to identify high versus low tumor GR expression.

Next we used the Kaplan–Meier method to estimate RFS for patients with NR3C1-high versus NR3C1-low expressing tumors (Fig. 3). We stratified patients into untreated and adjuvant-treated (tamoxifen treatment for ESR1⁺ patients and adjuvant chemotherapy for ESR1⁻ patients) groups. As hypothesized, we found that in ESR1⁻ tumors, high NR3C1 expression was indeed associated with significantly worse outcome in both the untreated (P = 0.001, log-rank test; HR = 2.23) and adjuvant chemotherapy–treated (P = 5.8e-7, log-rank test; HR = 6.83) groups. Because the analysis was not carried out prospectively, the predictive ability of GR expression (with respect to chemotherapy) in ER⁻ breast cancer cannot be assessed, but the data suggest that high GR expression is associated with a relatively poor prognosis regardless of adjuvant chemotherapy treatment (Fig. 3B and D). Unexpectedly, we observed the reverse association between high NR3C1 gene expression and clinical outcome in ESR1⁺ patients (Fig. 3A and C). In ESR1⁺ cancer, both untreated (P = 0.03, log-rank test; HR = 0.60) and tamoxifen-treated (P = 7.75e-8, log-rank test; HR = 0.25) patients with high NR3C1 tumors had better outcome compared with patients with low NR3C1-expressing tumors. Therefore, the phenotypic context of ER expression seems to reverse the association of high GR expression with a poor outcome in breast cancer.

To be certain that these results were not confounded by an association of NR3C1 expression with the expression of other receptors used in breast cancer prognosis, we also calculated the individual correlation coefficients between GR gene expression (216321_s_at) and genes encoding PR, androgen receptor (AR), and HER2. We found the correlation between GR expression and other receptors to be weak (Table 3), making it unlikely that expression of these receptors influenced our results. One of the HER2 probesets (210930_s_at) showed an approximately 30% negative correlation with NR3C1 expression. The influence of this group of HER2⁺ tumors on patient clinical outcome was therefore examined by excluding those patients in the top HER2 expression quartile from the ESR1⁻ dataset and then reanalyzing the remaining patients. Similar results, that is, a significantly worse prognosis in NR3C1-high tumors, were observed when we excluded the HER2-high patients (Supplementary Fig. S3), suggesting that HER2-positivity does not confound the association of high NR3C1 expression with a poor prognosis among ESR1⁻ patients.

We next examined the differences in global gene expression between high GR-expressing versus low GR-expressing tumors among ER⁺ and ER⁻ breast cancers as determined by tumor ESR1 expression. We carried out SAM (37) and identified genes showing a significant difference in expression (≥ 1.5-fold or ≤−1.5-fold change, false discovery rate ≤1%) in NR3C1-high versus NR3C1-low tumors (Fig. 4).

We found that in ESR1⁺ tumors, 4,069 genes were differentially expressed in NR3C1-high versus NR3C1-low tumors. Among these genes, 3,488 were significantly upregulated and 581 were downregulated. In ESR1⁻ tumors, 5,170 genes were significantly differentially expressed, among which 3,209 were upregulated and 1,961 were downregulated. We then classified these differentially regulated genes into 3 groups: (a) genes regulated exclusively in ESR1⁺ tumors (n = 930, 882 upregulated and 48 downregulated); (b) genes regulated exclusively in ESR1⁻ tumors (n = 2,031, 603 upregulated and 1,428 downregulated); and (c) shared genes commonly regulated in both ESR1 subtypes (n = 3,139, 2,606 upregulated and 533 downregulated; Fig. 4). Thus, the majority of genes associated with high NR3C1 expression were shared between ESR1⁻ and ESR1⁺ tumors, although there were important differences among the uniquely expressed genes. For example, the majority of the genes uniquely expressed in ESR1⁺ breast cancers were upregulated, whereas in ESR1⁻ tumors, the majority were downregulated, suggesting an influence of ER context on modulating the direction of GR-associated gene expression.

We then carried out a gene enrichment pathway analysis, using GeneGo's MetaCore software suite, to identify the most significant GR-associated pathways found in the ESR1⁺ and ESR1⁻ breast cancers (Fig. 4). Top pathways that were identified involved EMT, cell adhesion, and immune response. Somewhat surprisingly, we found that EMT-associated genes were differentially expressed in both the ESR1⁺ and ESR1⁻ breast cancer subtypes, suggesting that high GR expression influences EMT-associated gene expression in an ER-independent fashion. Many cell adhesion pathway genes were also identified in both ESR1⁺ and ESR1⁻ breast cancers, but the majority of genes in these pathways were upregulated in ESR1⁺ breast cancer and downregulated in ESR1⁻ tumors. Interestingly, a significant number of immune response pathway genes were only identified in ESR1⁻ breast cancers, suggesting that GR-associated modulation of the immune response may play a specific role in ER⁻ breast cancer outcome.

We next asked how the ER⁻ cell line data and ESR1⁻ human primary tumor data compared with respect to overlapping GR-associated gene expression. We examined genes that were common to both the 2,031 differentially regulated genes from ESR1⁻/GR-high tumors (Fig. 4) and the gene set (n = 187) derived from our cell line ChIP-seq and gene expression data. We found that 90 genes were identified in both analyses and likely represented direct GR target genes (Table 4). As hypothesized, these genes encoded antiapoptotic proteins previously identified as critical to GR-mediated tumor cell survival and chemotherapy resistance such as SGK1 and DUSP1/MKP1. Unexpectedly, genes encoding proteins central to EMT (e.g., SNAI2/SLUG), chromatin remodeling (e.g., SMARCA2), and epithelial cell–inflammatory cell interactions (e.g., IL1R1 and PTGDS) were also identified. Thus, the combined cell line and primary human tumor data analysis suggest that the GR directly mediates pathways associated with aggressive breast cancer biology. Taken together with the high recurrence rate found in ER⁻ patients with high GR-expressing early breast cancers, these finding suggest that GR expression may be a useful marker of high-risk early-stage ER⁻ breast cancer as well as a potential therapeutic target in these cancers.

In the clinical treatment of breast cancer, tumors are initially classified by ER, PR, and HER2 expressions. This classification guides both therapy and prognosis. Compared with early-stage ER⁺ breast cancer, similar stage ER⁻ breast cancer is more aggressive and relapses earlier. In addition, ER⁻ patients have relatively few options for adjuvant treatment because only ER⁺ patients benefit significantly from antiestrogen agents (5, 6). As a follow-up to our previous studies examining the molecular targets of GR activation that promote cell survival in ER⁻ breast epithelial and cancer cells (12–15), this study identified both genome-wide GR target genes in an ER⁻ breast epithelial cell line and GR-associated genes in early ER⁻ breast cancers. In our meta-analysis of 1,378 early-stage breast cancer patients, we found that patients with ER⁻/GR-high tumors have significantly increased risk of early relapse compared with patients with ER⁻/GR-low tumors. Unexpectedly, we also found that high GR expression was associated with better outcome in ER⁺ breast patients, suggesting potential cross-talk between the ER and GR (52, 53). Recently, loss of the E2F1-associated gene signature was also found to associate with a reversed prognosis in ER⁺ versus ER⁻ breast cancers, similarly suggesting the importance of ER expression in the activity and function of transcriptional networks in breast cancer (54). The unexpected finding of an improved RFS time for breast cancer patients with GR-high/ER⁺ tumors may in part be due to antagonism between ER and GR signaling (55). For example, in ER⁺ MCF-7 cells, glucocorticoid treatment inhibits the growth-stimulatory effect of estrogen (56, 57). Conversely, ER activation can also inhibit glucocorticoid action. Specifically, activated ER can induce the expression of the protein phosphatase 5 (PP5) gene, which in turn mediates the dephosphorylation and inactivation of the GR at Ser211 (52). Therefore, unlike ER⁻ tumors, where GR regulates genes independently of estrogen and ER-alpha action, in ER⁺ tumors, ER may alter GR action and the GR may in turn also alter a subset of important ER-regulated genes involved in tumor growth. This hypothesis will be explored in further studies examining the role of ligand-activated ER-alpha expression in GR chromatin occupancy and transcriptional activity.

Analysis of gene expression pathways, based on our meta-dataset of human tumors, revealed that the EMT pathway is associated with high GR expression in both ER⁻ and ER⁺ tumors, suggesting that GR expression may mediate transcriptional activation of fundamental cytoskeletal genes. Interestingly, however, cell adhesion pathway gene expression was differently regulated between ER⁻ and ER⁺ tumor subtypes; in ER⁺ tumors, genes encoding cell adhesion proteins were generally upregulated whereas in ER⁻ tumors they were downregulated. In contrast, inflammatory and macrophage-associated genes (i.e., immune response) were uniquely identified as differentially expressed only in the ER⁻ tumors. Furthermore, the upregulation of several genes encoding proteins involved in epithelial cell–macrophage attraction such as DAP12 and Syk was observed in the gene set common to both the ER⁻ cell line data and ER⁻ primary tumor gene expression. This interesting finding suggests that GR activity plays an important role in tumor epithelial cell communication with the microenvironment in the ER⁻ breast cancer subtype (58).

On the basis of these findings, the role of GR overexpression in ER⁻ breast cancer should be evaluated prospectively in trials of early-stage patients undergoing adjuvant therapy. Eventually, the presence of high GR expression might be a useful tool to determine a recurrence threshold for employing systemic treatment in ER⁻ patients with less than 1.0 cm node-negative breast cancer who currently do not routinely receive adjuvant systemic therapy but appear likely to relapse if they have the ER⁻/GR⁺ subtype. Furthermore, GR expression may represent a novel therapeutic target for chemotherapy-resistant ER⁻ breast cancer that should be explored with the new second generation specific GR antagonists currently under development.

No potential conflicts of interest were disclosed.

We thank Drs. Gini Fleming and Geoffrey Greene for valuable discussions and Dr. John Foekens for sharing unpublished data.

This study was supported by the AVON Foundation for Women, NIH R01CA089208 (S.D. Conzen), pilot project funding from P50CA125183 (SPORE), and the University of Chicago Comprehensive Cancer Center Core Facility Support P30CA014599-35.

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.