This article is mentioned in:

Introduction

Breast cancer is a common types of malignancy among females, worldwide (1). Genetic and epigenetic alterations are involved in the underlying mechanisms associated with breast cancer development (2–4). Although therapeutic and diagnostic methods have improved, this type of cancer remains the primary cause of cancer-associated mortality among females (5). It is estimated that there are 464,000 cases of breast cancer, accounting for 13.5% of all cancer cases in Europe in 2012, and the number of breast cancer-associated mortalities is 131,000 (6). Furthermore, breast cancer is the most common cause of cancer-associated mortality in females. Therefore, it is essential to understand its molecular mechanism and develop more effective therapeutic methods for breast cancer treatment.

The estrogen receptor (ER) is critical in determining the phenotype of human breast cancers and is one of the most important therapeutic targets (7). Furthermore, certain studies have suggested that activation of ER is responsible for various biological processes, including cell growth and differentiation, and programmed cell death (8,9). It is reported that the response of ER to estrogen is critical in controlling specific protein synthesis (10). ER-mediated transcription has been extensively investigated on a small number of endogenous target promoters (11,12). Carroll et al (13) identified various novel features of ER transcription, including an involvement of distal cis-regulatory enhancer regions, and a requirement for the Forkhead protein, FoxA1, in facilitating ER binding to chromatin and subsequent gene transcription. However, the mechanisms underlying estrogen-associated gene expression changes in breast cancer remain poorly understood.

The human cancer cell line, MCF7, contains ER and demonstrates an estrogen response (14). In the present study, the microarray data was obtained from the Gene Expression Omnibus (GEO) database, which was developed using MCF7 cells stimulated with estrogen for different durations. The differentially expressed genes (DEGs) between the control and estrogen treatment groups were analyzed. The ER binding sites were identified and motif analysis was performed. The aim was to investigate the regulatory mechanism of ER in the progression of breast cancer.

Materials and methods

Affymetrix microarray data

The mRNA microarray datasets (accession no. GSE11324) were obtained from the GEO database (http://www.ncbi.nlm.nih.gov/geo/), which was deposited by Carroll et al (13). The gene expression profiles were developed from 12 batches of MCF7 cells that had been treated with 100 nM estrogen for 0, 3, 6 or 12 h. The experiment had been repeated three times.

Data preprocessing and DEG Identification

The raw Affymetrix CEL data were downloaded based on the platform of GPL570 (HG-U133_Plus_2) Affymetrix GeneChip® Human Genome U133 Plus 2.0 Array (Affymetrix, Santa Clara, CA, USA). To obtain the expression matrix, mRNA expression data were first preprocessed using the Affy Package in R language (http://bioconductor.org/packages/release/bioc/html/affy.html) (15). When multiple probes mapped to the same Entrez Gene ID, the mean expression value of these probes was calculated for the gene.

The DEGs between the ER-stimulated MCF7 cells were screened by limma package in R language (16) and genes with an adjusted P-value <0.05 were considered to be DEGs. Subsequently, the DEGs in the 3 vs. 0 h, 6 vs. 0 h and 12 vs. 0 h groups underwent hierarchical clustering analysis (17) using the pheatmap package in R language (R Core Team, Vienna, Austria).

ER binding site analysis

Chromatin immunoprecipitation with high-throughput sequencing (ChIP-seq) data was obtained from the GEO repository (accession no. GSE25710; ChIP-seq for forkhead box A1, ER and CCCTC-binding factor in breast cancer cell lines). The ChIP-seq data were treated with ER antibody (Ab-10; Neomarkers, Lab Vision; Fremont, CA, USA) and mapped to the whole genome sequence using BowTie software version 2.1.0 (https://sourceforge.net/projects/bowtie-bio/files/bowtie2/) (18). To identify the possible ER binding sites, peak calling was performed by model-based analysis of ChIP-Seq (MACS; version 1.4.2) (19). The q-value was set at <0.01 and served as the cut-off to improve the ChIP-seq peak detection. The DNA sequence fragments that interacted with ER in different regions of the whole genome were evaluated using the Cis-regulatory Element Annotation System software (http://liulab.dfci.harvard.edu/CEAS) and a graph of the results was constructed (20,21).

Binding and expression target analysis (BETA)

BETA version 1.0.7 (http://cistrome.org/BETA/) (22) is a software package that predicts target genes through the analysis of ChIP-seq data and DEGs. The ER target genes were screened, using BETA software, based on the ER-associated ChIP-seq data and ER-stimulated DEG data. Briefly, the genes within a 100-kb distance of the significant top 10,000 peaks were collected. These genes were ranked based on their distance from the peak and the significance of their differential expression. The genes ranked at the top exhibited the greatest possibility of being regulated by ER.

Subsequently, the regulatory function (activation or inhibition) of ER among the top 500 DEGs in the 12 vs. 0 h group, which neared the top 10,000 peaks were evaluated by BETA analysis. The motifs of the ER binding sites were obtained and the factors that interacted with ER were predicted (23).

Results

Identification of DEGs

The DEGs with adjusted P<0.05 in the 3 vs. 0 h, 6 vs. 0 h and 12 vs. 0 h groups were screened out. Subsequently, the gene expression profiles of DEGs at different time-points were analyzed. Hierarchical clustering indicated that the DEGs were clearly separated and the difference in gene expression became more pronounced with increasing duration of estrogen treatment (Fig. 1). In order to obtain more reliable results, a total of 3,122 DEGs with adjusted P<0.01 in the 12 vs. 0 h group were selected for further analysis, including 1,755 upregulated and 1,366 downregulated genes.

Figure 1 Heatmap showing the results of hierarchical cluster analysis of DEGs at four different time-points. The gene expression levels from high to low are indicated in red to green, respectively, and the number of DEGs is expressed in pink to blue as indicated. DEG, differentially expressed gene.

Identification of ER-specific binding sites

Based on the ChIP-seq peaks determined by MACS and the cut-off value of q<0.01, a total of 10,058 peaks were obtained. Based on the ChIP-seq data analysis, the distribution of ER binding sites in the whole genome was determined. Fig. 2 demonstrates that the DNA sequence fragments that interact with the ER are located in different regions of the whole genome, such as promoter, downstream and intergenic regions.

Figure 2 Distribution of ER binding sites in the genome. (A) Prediction of ER binding sites in the whole genome. (B) ER binding sites tested by chromatin immunoprecipitation sequencing, indicating a higher number of ER binding sites in the promoter region compared with A. ER, estrogen receptor; UTR, untranslated region.

BETA analysis

In order to analyze the regulatory effects of ER on its target genes, the 10,000 most significant peaks and the 500 most significant DEGs in the 12 vs. 0 h group were selected for BETA analysis. Results of the BETA analysis indicated that the ER exerts an inhibitory role in addition to an activating role regarding the regulation of its target genes (Fig. 3). In addition, motif analysis revealed that ER may exert an activation role in gene expression by interacting with MEIS1 and FOXP3 (Table I, Part A) and inhibit gene expression by interacting with THRB and GRHL1 (Table I, Part B).

Figure 3 ER activating/repressive function prediction. Genes were cumulated by ranking on the basis of the regulatory potential score from high to low. A significance of ER activation/inhibition in gene expression was indicated. ER, estrogen receptor.

Table I Motifs in the target genes.

Table I

Motifs in the target genes.

A, Upregulated target genes
Symbol	DNA-binding domain	Species	P-value (t-test)	T-score
ESR1	Nuclear hormone receptor family	Homo sapiens	2.24×10−40	14.50
ESRRB	Nuclear hormone receptor family	Homo sapiens	2.24×10−40	14.50
ESRRA	Nuclear hormone receptor family	Homo sapiens	2.24×10−40	14.50
ESRRG	Nuclear hormone receptor family	Homo sapiens	2.24×10−40	14.50
RARA	Nuclear hormone receptor family	Homo sapiens	2.24×10−40	14.50
PPARG	Nuclear hormone receptor family	Homo sapiens	2.24×10−40	14.50
NR2F1	Nuclear hormone receptor family	Homo sapiens	2.24×10−40	14.50
ESR2	Nuclear hormone receptor family	Homo sapiens	2.24×10−40	14.50
MEIS1	Homeo domain family	Homo sapiens	1.90×10−18	8.95
FOXP3	Forkhead domain family	Homo sapiens	2.52×10−9	5.93

B, Downregulated target genes
Symbol	DNA-binding domain	Species	P-value (t-test)	T-score
ESR1	Nuclear hormone receptor family	Homo sapiens	1.12×10−43	14.58
ESRRB	Nuclear hormone receptor family	Homo sapiens	1.12×10−43	14.58
ESRRA	Nuclear hormone receptor family	Homo sapiens	1.12×10−43	14.58
ESRRG	Nuclear hormone receptor family	Homo sapiens	1.12×10−43	14.58
RARA	Nuclear hormone receptor family	Homo sapiens	1.12×10−43	14.58
PPARG	Nuclear hormone receptor family	Homo sapiens	1.12×10−43	14.58
NR2F1	Nuclear hormone receptor family	Homo sapiens	1.12×10−43	14.58
ESR2	Nuclear hormone receptor family	Homo sapiens	1.12×10−43	14.58
THRB	Nuclear hormone receptor family	Homo sapiens	7.83×10−12	6.82
GRHL1	CP2 transcription factor domain family	Homo sapiens	1.28×10−11	6.75

[i] ESR1, estrogen receptor 1; ESRRB, estrogen related receptor β; ESRRA, estrogen related receptor α; ESRRG, estrogen related receptor γ; RARA, retinoic acid receptor α; PPARG, peroxisome proliferator activated receptor γ; NR2F1, nuclear receptor subfamily 2 group F member 1; ESR2, estrogen receptor 2 (ER β); MEIS1, meis homeobox 1; FOXP3, forkhead box P3; THRB, thyroid hormone receptor, β; GRHL1, grainyhead-like 1.

Discussion

The ER is recognized as the master transcriptional regulator of the breast cancer phenotype and is critical in predicting the early recurrence of breast cancer. However, to the best of our knowledge, the role of ER in breast cancer cell gene expression has not been clearly clarified. In the present study, DEGs in the MCF7 breast cancer cell line that had been stimulated by estrogen for different durations were analyzed. Using a combination of the ChIP-seq dataset and the identified DEGs, the ER target and response genes were predicted. A set of cis-acting targets across the whole genome of the ER were identified.

The present results demonstrated that 3,122 genes were differentially expressed as a result of estrogen stimulation. The hierarchical clustering analysis for the DEGs indicated that the long-term stimulation by estrogen improved the differential gene expression in breast cancer cells. Using motif analysis, ER was identified to inhibit and stimulate target gene expression, which was demonstrated by interactions between MEIS1, FOXP3, THRB and GRHL1, and ER.

In the present study, ER was found to activate gene expression by interacting with MEIS1 and FOXP3. MEIS1 is a homeobox gene, and encodes the homeobox protein, MEIS1 (24). In addition, certain homeobox proteins are associated with tumor formation; Mahmoud et al (25) proposed MEIS1 as a critical transcriptional regulator of cardiomyocyte proliferation and as a potential therapeutic target for heart regeneration. It was previously reported that MEIS1 is a prognostic and predictive biomarker for breast cancer (26) and a recent study indicated that MEIS1, as a HOX gene, is associated with decreased proliferation in the mesenchymal stem-like subtype of breast cancer (27). Furthermore, the androgen/estrogen metabolism pathway is responsible for ER-negative breast cancer (27). Therefore, whether the MEIS1 response to ER is responsible for the ER-positive breast cancer subtype requires further investigation.

FOXP3, a member of the FOX protein family, is involved in the immune system response. FOXP3 controls the expression of numerous genes and has recently been reported to be expressed in tumor cells (28). FOXP3 expression was reported to be enhanced in estrogen-treated mice (29). Fox et al (30) showed that expression of the transcription factor, FOXP1 is associated with ERα and improved survival in patients with primary breast carcinomas. Merlo et al (28) suggested that FOXP3 expression was a novel independent prognostic factor for breast carcinoma. Thus, ER may interact with MEIS1 and FOXP3 to activate gene expression in breast cancer. In addition, the results of the present study showed that ER suppressed gene expression via THRB and GRHL1. THRB is considered to be a potential cancer suppressor (31) and THRB gene silencing by aberrant methylation is highly prevalent in breast cancer patients (31). Furthermore, the loss of THRB expression as a result of methylation may be a plasma biomarker for the prognosis of breast cancer patients (31). Baniahmad et al (32) identified that the interaction of THRB with transcription factors may mediate the activation of the target gene. GRHL1 inhibits tumorigenicity and is a prognostic marker in neuroblastoma (33,34). Tao et al (35) found that Xenopus GRHL1 is essential for epidermal differentiation. de la Garza et al (36) indicated that interferon regulatory factor 6 promoted differentiation of the periderm by stimulating the expression of grainyhead-like 3 (36). Although the evidence for the interaction between ER and THRB and GRHL1 is insufficient, THRB and GRHL may be potential targets to analyze the function of ER in breast cancer.

In conclusion, ER may activate or suppress gene expression by interacting with MEIS1 and FOXP3, or THRB and GRHL1, respectively. These data may be useful for identifying novel therapeutic agents and designing clinical trials; however, further experiments are required to confirm the results.

References