Profiling of gene expression regulated by 17β-estradiol and tamoxifen in estrogen receptor-positive and estrogen receptor-negative human breast cancer cell lines

Abstract

One area of great importance in breast cancer (BC) research is the study of gene expression regulated by both estrogenic and antiestrogenic agents. Although many studies have been performed in this area, most of them have only addressed the effects of 17β-estradiol (E2) and tamoxifen (TAM) on MCF7 cells. This study aimed to determine the effect of low doses of E2 and TAM on the expression levels of 84 key genes, which are commonly involved in breast carcinogenesis, in four BC cell lines differentially expressing estrogen receptor (ER) α and HER2 (MCF7, T47D, BT474, and SKBR3). The results allowed us to determine the expression patterns modulated by E2 and TAM in ERα+ and ERα− cell lines, as well as to identify differences in expression patterns. Although the MCF7 cell line is the most frequently used model to determine gene expression profiles in response to E2 and TAM, the changes in gene expression patterns identified in ERα+ and ERα− cell lines could reflect distinctive properties of these cells. Our results could provide important markers to be validated in BC patient samples, and subsequently used for predicting the outcome in ERα+ and ERα− tumors after TAM or hormonal therapy. Considering that BC is a molecularly heterogeneous disease, it is important to understand how well, and which cell lines, best model that diversity.

Keywords:

• breast cancer
• cell lines
• 17β-estradiol
• tamoxifen
• ERα+
• ERα−
• qPCR

Introduction

Breast cancer (BC) is the most common cancer in women,¹ which has been associated with high exposure to 17β-estradiol (E2). E2 is the main estrogenic hormone, playing an important regulatory role in biological processes such as mammary epithelial cell proliferation and breast development during sexual maturity.² The pleiotropic effects of E2 on its target tissues are mediated via its receptors, estrogen receptor (ER) α and β, which are members of the nuclear receptor superfamily of ligand-activated transcription factors.³ E2 induces well-established mechanisms to either activate or suppress transcription of its target genes, such as direct interaction with DNA at estrogen response elements,⁴ through transcription factor cross talk,⁵ and non-genomic mechanisms whereby E2 binds to ERs localized in the cell membrane and activates signal transduction pathways.⁶

The ability of ERs to regulate transcription is also dependent on the nature of their ligands with various natural and synthetic selective ER modulators acting as either ER agonists or ER antagonists.^7,8 One of them is tamoxifen (TAM), a well-known nonsteroidal antiestrogen with partial agonistic activity, which is extensively used in the treatment of ERα+ BC. Although ERα positivity is a well-established predictor of responses to TAM with ERα− patients considered as nonresponders, it is known that 5–10% of ERα− tumors benefit from adjuvant TAM treatment.^9–11 However, responses to TAM have frequently limited durations due to the development of resistance.^12,13

The MCF7 BC cell line is one of the most widely used models to study estrogen signaling, because the cells express high levels of ERα and their growth is strongly stimulated by E2 in vitro and in xenografts models.^14,15 Considering these aspects, several studies have evaluated gene expression profiles, which focused on identifying changes in gene expression associated with BC prognoses,^16–19 drug resistance,²⁰ or tumor aggressiveness.²¹ Nonetheless, relatively few studies have evaluated gene expression profiles in other ERα+ (T47D and BT474) or ERα− (SKBR3) BC cell lines. Similarly, the number of studies that have evaluated the effect of TAM on the gene expression profiles in these cell lines is limited.

Considering that different cell lines with similar immunophenotypes are used in basic and clinical BC research, identification and comparison of genes and genetic pathways responsive to E2 and TAM are necessary to understand the complex effects of these agents on BC cells. This study aimed to determine the effects of low doses of E2 and TAM (0.01 and 1 µM, respectively) on the expression levels of 84 key genes, which are commonly involved in the regulation of signal transduction and other biological processes related to breast carcinogenesis in ERα+ (MCF7, T47D, and BT474) and ERα− (SKBR3) BC cell lines.

The expression patterns of E2- and TAM-regulated genes identified in ERα+ and ERα− cell lines allowed the identification of variations in gene responses to E2 and suggested genes that might serve as markers of sensitivity and/or resistance to TAM. We found that the expression patterns of some of the genes evaluated in this study clearly discriminate ERα+ and ERα− BC cell lines.

Materials and methods

Cell culture

Human BC cell lines MCF7 (ERα+/HER2−), T47D (ERα+/HER2−), BT474 (ERα+/HER2+), and SKBR3 (ERα−/HER2+) were obtained from the American Type Culture Collection (Manassas, VA, USA). MCF7, T47D, and SKBR3 cells were cultured in RPMI 1640 medium (Sigma, St. Louis, MO, USA), whereas BT474 cells were cultured in Dulbecco’s modified Eagle’s medium (Sigma). All culture media were supplemented with 10% fetal bovine serum (FBS) (Sigma), antibiotic–antimycotic solution (1X) (Sigma), and 2 mM l-glutamine (Thermo Fisher Scientific, Waltham, MA, USA). Cells were maintained at 37°C and 5% CO₂. Absence of mycoplasma contamination was verified by conventional PCR assays.

E2 and TAM treatments

At 48 h prior to the addition of E2 (E2758; Sigma) or TAM (T5648; Sigma), culture media were changed to phenol red-free RPMI 1640 (Sigma) containing 10% charcoal-stripped FBS (Sigma) to remove endogenous serum steroids and eliminate the weak estrogen agonistic activity of phenol red.²² E2 and TAM were dissolved in absolute ethanol, diluted in the medium at 0.01 and 1 µM, respectively, and then applied to cells for 24 and 48 h. Each time course experiment was carried out in triplicate. RNA was extracted in all cases and used for cDNA synthesis and quantitative real-time PCR (qPCR) assays. Cells without treatment were used as controls.

RNA extraction

Cell lines were collected in 1 ml TRIZOL Reagent^® (Ambion; Thermo Fisher Scientific), and RNA extraction was performed according to the manufacturer’s instructions. Briefly, total RNA was separated by addition of chloroform, followed by precipitation with isopropanol and 75% ethanol. RNA pellets were resuspended in diethyl pyrocarbonate-treated water (Ambion; Thermo Fisher Scientific) and quantified by spectrophotometry using a NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific). RNA was further purified using RNeasy MinElute Cleanup columns (SABiosciences™; Qiagen, Hilden, Germany) and treated with an RNase-Free DNase Set (SABiosciences™; Qiagen) to remove contaminating genomic DNA. Purified RNA samples were stored at −80°C until analysis.

Reverse transcription and qPCR

One microgram of RNA, extracted from cell lines treated with either E2 or TAM for 24 and 48 h, was used to synthesize cDNA with an RT² First Strand Kit (SABiosciences™; Qiagen) according to the manufacturer’s instructions. The cDNA for each qPCR array test was combined with ready-to-use RT² SYBR Green qPCR Master Mix (SABiosciences™; Qiagen), and equal aliquots of this mixture (25 µL) were added to each well of a Human Breast Cancer RT² Profiler™ PCR Array plate (PAHS-131A, SABiosciences™; Qiagen). This PCR array plate contained gene-specific primer sets for 96 genes including 84 genes associated with breast carcinogenesis, five reference genes (B2M, HPRT1, RPL13A, GAPDH, and ACTB), one genomic DNA control, three reverse transcription controls, and three positive PCR controls (Table 1). cDNA was amplified using an iCycler^® Real-Time PCR System (Bio-Rad Laboratories, Hercules, CA, USA) for 40 cycles (95°C for 15 s and 60°C for 1 min). The threshold cycle (Ct) values were used to calculate the fold changes (FCs) in gene expression (treated cells vs controls – untreated cells) using the PCR array data analysis software (SABiosciences™; Qiagen), which was based on the 2^−ΔΔCt method. Gene expression was classified as significantly regulated by E2 or TAM treatments when the FC in expression was at least ±2. qPCR experiments were carried out in three biological replicates, and from each biological replicate, two technical replicates were also performed. The results shown are the average of all the experiments performed.

Table 1 Array layout of genes in the Human Breast Cancer RT² Profiler™ PCR Array plate (PAHS-131A, SABiosciences™)

ABCB1	ABCG2	ADAM23	AKT1	APC	AR	ATM	BAD	BCL2	BIRC5	BRCA1	BRCA2
CCNA1	CCND1	CCND2	CCNE1	CDH1	CDH13	CDK2	CDKN1A	CDKN1C	CDKN2A	CSF1	CST6
CTNNB1	CTSD	EGF	EGFR	ERBB2	ESR1	ESR2	FOXA1	GATA3	GLI1	GRB7	GSTP1
HIC1	ID1	IGF1	IGF1R	IGFBP3	IL6	JUN	KRT18	KRT19	KRT5	KRT8	MAPK1
MAPK3	MAPK8	MGMT	MKI67	MLH1	MMP2	MMP9	MUC1	MYC	NME1	NOTCH1	NR3C1
PGR	PLAU	PRDM2	PTEN	PTGS2	PYCARD	RARB	RASSF1	RB1	SERPINE1	SFN	SFRP1
SLC39A6	SLIT2	SNAI2	SRC	TFF3	TGFB1	THBS1	TP53	TP73	TWIST1	VEGFA	XBP1
B2M	HPRT1	RPL13A	GAPDH	ACTB	HGDC	RTC	RTC	RTC	PPC	PPC	PPC

Heatmaps of differentially expressed genes (DEGs)

Significant (P<0.05) DEGs (FC: ±2) were selected to graphically evaluate the activation or inhibition of gene expression in each treated cell line. Genes were clustering using the gplots package from the Bioconductor project.²³ Euclidean distance was used to calculate the matrix of distances, and clusters were generated using Ward’s method.

Biological interpretation of the DEGs (pathway analysis)

Pathway enrichment analysis was conducted by over-representation analysis (ORA). ORA was performed in DAVID Bioinformatics Resources 6.7 of the National Institute of Allergy and Infectious Diseases (National Institutes of Health, Rockville, MD, USA)²⁴ by employing a pre-filtered list of DEGs. Overrepresented pathways were generated based on the information in the Gene Ontology (GO) project compared with a background limited to the genes analyzed in this study (84 genes). Fisher’s exact test was performed to determine the likelihood of obtaining at least equivalent numbers of genes by an actual overlap between the input gene set and the genes in each identified GO category.

Overall representation of the data

To obtain a simultaneously general representation of the gene expression response of the variously treated cell lines, principal component analysis (PCA) was performed with FCs of all evaluated genes. Because the FC is a ratio of the expression level with a specific treatment and the expression level without treatment, a reference time point for all experimental groups was created (control time 0, FC=1 for all genes). PCA is a linear projection method and an unsupervised exploratory technique used to remove noise, reduce dimensionality, and identify common/dominant signals oriented to find biological meaning. The two principal components (PCs) with the highest degree of variance were plotted. PCA was performed with the prcomp package, and the plot was drawn with gplots, both from the Bioconductor project.²³

To identify the main genes responsible for the different PC segments, we used the pcaGoPromoter package²⁵ implemented in R and part of the Bioconductor project. This methodology allows identification of the loading genes in each PC. Additionally, it quantifies the importance of the genes in the reduced predictor. Hansen et al²⁵ empirically showed that between 0.5 and 8.5% of the total genes is the window to perform the selection of loading genes for good and stable ORA. In this study, a 5% threshold was selected (five genes).

Statistical analysis

Student’s t-test was performed to determine significant differences in FC values between groups (treated and untreated cells). P-values of <0.05 (*) and 0.01 (**) were considered as statistically significant. All tables and graphs show the mean values of three independent experiments.

Results

Gene expression patterns in BC cell lines treated with E2

In MCF7 cells, 19 out of 84 analyzed genes were regulated by E2 at 24 and/or 48 h (22.6%). Among these genes, 47.4% (n=9) were induced and 52.6% (n=10) were suppressed (Table 2 and Figure 1A). These modulated genes were assigned to one of three patterns: genes regulated at 24 h only (early expression; 15.8%), genes regulated at 24 and 48 h (early and late expression; 68.4%), and genes regulated at 48 h only (late expression; 15.8%). Cluster analysis demonstrated two patterns of modulated gene expression at both time points: upregulated (first cluster) and downregulated (second cluster) (Figure 1A). Among the 19 genes regulated by E2 in MCF7 cells, BIRC5, NME1, BRCA1, BRCA2, EGF, HPRT1, MKI67, PGR, and TFF3 were upregulated significantly (P<0.05), while CCND2, IGFBP3, RARB, SLIT2, NR3C1, and CDKN1A genes were downregulated significantly (Table 2).

Figure 1 Cluster analysis of the time course of E2-regulated gene expression in (A) MCF7, (B) T47D, (C) BT474, and (D) SKBR3 cells. Gene cluster analysis was performed for 84 genes after E2 exposure at 24 and 48 h. The threshold cycle (Ct) values were used to calculate the FCs in gene expression (treated cells vs controls – untreated cells) using the PCR array data analysis software (SABiosciences™), which was based on the 2^−ΔΔCt method. Significant differentially expressed genes (FC=±2; P<0.05) were selected to graphically evaluate the activation or inhibition of gene expression in each treated cell line. Induced genes are shown in red, suppressed genes in green, and unregulated genes in black.

Abbreviations: E2, 17β-estradiol; FC, fold change.

Table 2 List of E2-regulated genes in MCF7, T47D, BT474, and SKBR3 cells

Cell line	24 h	48 h	Gene name	GenBank no.	Pattern expression
MCF7	3.4822^**	1.8025	BIRC5	NM_001168	Early
	2.1435^*	1.7411	NME1	NM_000269
	−2.4623	−1.2311	SERPINE1	NM_000602
	2.4623^**	2.639^**	BRCA1	NM_007294	Early and late
	3.249^**	2.3784^**	BRCA2	NM_000059
	−2^**	−2.0705^**	CCND2	NM_001759
	3.4822^*	5.0982^**	EGF	NM_001963
	2.8284^**	2.4623^**	HPRT1	NM_000194
	−3.0314	−2.0705	TWIST1	NM_000474
	−2.4623	−2.4623	ERBB2	NM_004448
	−2.8945	−3.1383	GRB7	NM_005310
	−2.639^**	−5.0982^**	IGFBP3	NM_000598
	3.249^**	2.0705^**	MKI67	NM_002417
	6.0629^**	8.8766^**	PGR	NM_000926
	−2.1435^**	−2.2191^**	RARB	NM_000965
	−2.639^**	−3.7321^**	SLIT2	NM_004787
	1.7613	2.5491^**	TFF3	NM_003226	Late
	−1.7211	−2.1435^**	NR3C1	NM_000176
	−1.7211	−2.2974^**	CDKN1A	NM_000389
T47D	2.9794^**	1.9656	BRCA1	NM_007294	Early
	2.7678^*	1.9571	CCNE1	NM_001238
	2.1067^**	1.7715	PGR	NM_000926
	2.9409	1.8742	TP73	NM_005427
	2.6735	2.3274	ABCG2	NM_004827	Early and late
	2.0885^*	2.4837^*	BIRC5	NM_001168
	2.3784^**	5.278^**	MYC	NM_002467
	2.2482^**	2.4095^**	BRCA2	NM_000059
	2.3784	2.3784	SFN	NM_006142
	2.2875	2.2095	SRC	NM_005417
	2.7203	3.716^**	ESR2	NM_001437
	1.858	2.7203^**	MKI67	NM_002417	Late
	1.5223	2.0795	NOTCH1	NM_017617
	1.9487	2.0885	RASSF1	NM_007182
	1.3779	3.3782^*	TGFB1	NM_000660
	1.3899	2.2579^**	THBS1	NM_003246
	1.5422	2.2776	CTSD	NM_001909
BT474	−2.5669^*	−1.3013	SNAI2	NM_003068	Early
	2.2815^**	1.8277	NME1	NM_000269
	3.1167^**	3.6553^**	THBS1	NM_003246	Early and late
	−2.5669	−2.7895	GATA3	NM_002051
	3.1167^**	2.8679^**	IGF1R	NM_000875
	−2.4794^**	−3.9449^**	EGFR	NM_005228
	−4.9588	−2.8879	IGFBP3	NM_000598
	25.8125^**	15.1369^**	MYC	NM_002467
	9.1261^**	9.3179^**	PGR	NM_000926
	−2.3134^*	−2.6945^*	MUC1	NM_001018016
	−2.2346	−2.042	SLIT2	NM_004787
	−1.7532	−2.514	ERBB2	NM_004448	Late
	−1.3287	−2.114	ESR1	NM_000125
	−1.4743	−2.2658	FOXA1	NM_004496
	−1.5263	−2.114^*	GRB7	NM_005310
	−1.0792	−2.042^*	GSTP1	NM_000852
	−1.815	−2.042^**	NR3C1	NM_000176
	1.3566	4.0558^**	ABCG2	NM_004827
SKBR3	2.0849	1.1096	BAD	NM_004322	Early
	2.0849^*	1.1096	BRCA2	NM_000059
	3.7581^**	1.366	CCNE1	NM_001238
	3.2716	1.366	IGFBP3	NM_000598
	5.1337^**	1.4142	IL6	NM_000600
	−2.4453^**	−1.2746	MUC1	NM_001018016
	2.8481^*	1.9319	PLAU	NM_002658
	−2.1287	−1.0718	PYCARD	NM_013258
	2.395	1.1096	SERPINE1	NM_000602
	−2.1287^*	−2.2974^**	CDKN1C	NM_000076	Early and late
	4.0278	2.1435^**	BCL2	NM_000633
	1.3566	2.2191^**	CDH1	NM_004360	Late

Note: Values are fold changes of ±2 at each time point relative to control with the statistically significant values indicated by asterisks (^*P<0.05 and ^**P<0.01).

Abbreviation: E2, 17β-estradiol.

In T47D cells, expressions of 17 out of the 84 analyzed genes were modulated at 24 and/or 48 h (20.2%). In contrast to MCF7 cells, all E2-regulated genes were upregulated (Table 2 and Figure 1B). Among these genes, four (23.5%) showed early expression, six (35.3%) showed early and late expression, and seven (41.25%) showed late expression. Cluster analysis demonstrated three patterns of modulated gene expression with the first cluster including genes with early and late expression, the second cluster including genes regulated at both 24 and 48 h, and the third cluster corresponding to genes mostly regulated at 48 h (Figure 1B). Significantly altered expression (P<0.05) was observed in 10 (BRCA1, CCNE1, PGR, BIRC5, MYC, BRCA2, ESR2, MKI67, TGFB1, and THBS1) out of 17 genes (58.8%) upregulated by E2 in T47D cells (Table 2).

In BT474 cells, 18 out of the 84 analyzed genes (21.4%) were regulated by E2 at 24 and/or 48 h. Six (33.3%) were induced, while 12 (66.6%) were suppressed (Table 2 and Figure 1C). Among these genes, two (11.11%) showed early expression, nine (50%) showed early and late expression, and seven (38.88%) showed late expression. Similar to MCF7 cells, cluster analysis demonstrated two patterns of modulated gene expression with the first cluster including upregulated genes and the second cluster including downregulated genes at both examined time points (Figure 1C). However, it is important to note that, in contrast to MCF7 cells, most genes (six out of seven) were suppressed at 48 h (late expression). Twelve genes were regulated significantly (P<0.05) by E2 in BT474 cells: NME1, THBS1, IGF1R, MYC, PGR, and ABCG2 were upregulated, while SNAI2, EGFR, MUC1, GRB7, GSTP1, and NR3C1 were downregulated (Table 2).

In SKBR3 (ERα−) cells, E2 treatment resulted in the lowest number of modulated genes, 12 out of 84 (14.3%). Among them, nine (75%) were induced and three (25%) were suppressed (Table 2 and Figure 1D). In contrast to the other cell lines, most genes (nine) underwent early regulation at 24 h only. Cluster analysis demonstrated three patterns of modulated gene expression: upregulated genes (FC: >3) with an early response, upregulated genes (FC: <3) with an early response, and downregulated genes with an early response (Figure 1D). Among the 12 genes regulated by E2, six were significantly upregulated (BRCA2, CCNE1, IL6, PLAU, BCL2, and CDH1), while the significantly downregulated genes included MUC1 and CDKN1C (Table 2).

Gene expression patterns in BC cell lines treated with TAM

The number of TAM-regulated genes was lower compared with the E2 response in all cell lines. In MCF7 cells, five out of 84 analyzed genes (5.95%) changed their pattern of expression at 24 and/or 48 h: two of them were induced, while 3 were suppressed (Table 3 and Figure 2A). Only in MCF7 cells, at least one gene was seen in each of the three patterns: three genes showed early expression, one gene showed early and late expression, and one gene showed late expression. Only PGR and EGF (40%) showed significant increases in expression (P<0.05).

Figure 2 Cluster analysis of the time course of TAM-regulated gene expression in (A) MCF7, (B) T47D, (C) BT474, and (D) SKBR3 cells. Gene cluster analysis was performed for 84 genes after TAM exposure at 24 and 48 h. The threshold cycle (Ct) values were used to calculate the FCs in gene expression (treated cells vs controls – untreated cells) using the PCR array data analysis software (SABiosciences™), which was based on the 2^−ΔΔCt method. Significant differentially expressed genes (FC=±2; P<0.05) were selected to graphically evaluate the activation or inhibition of gene expression in each treated cell line. Induced genes are shown in red, suppressed genes in green, and unregulated genes in black.

Abbreviations: TAM, tamoxifen; FC, fold change.

Table 3 List of genes up- or downregulated (negative values) by TAM in MCF7, T47D, BT474, and SKBR3 cells

Cell line	24 h	48 h	Gene name	GenBank no.	Pattern expression
MCF7	−2.3784	−1.7613	CDKN1C	NM_000076	Early
	−2.0705	−1.6434	CSF1	NM_000757
	3.1383^*	1.977	PGR	NM_000926
	4.2871^**	3.6893^**	EGF	NM_001963	Early and late
	−1.366	−2.0232	SERPINE1	NM_000602	Late
T47D	2.3681	1.9235	GATA3	NM_002051	Early
	2.4516	1.3601	RARB	NM_000965
	2.2677	1.2152	TWIST1	NM_000474
BT474	−1.2746	−2.4284^**	EGFR	NM_005228	Late
	−1.5157	−2.6027	ERBB2	NM_004448
	−1.7411	−2.3457	GATA3	NM_002051
	−1.1892	−2.114	GRB7	NM_005310
	−1.1096	−2.114^**	JUN	NM_002228
	−1.2311	−2.3457	NOTCH1	NM_017617
	1.8025	2.6759^*	PGR	NM_000926
	−1.4641	−2.6945	RASSF1	NM_007182
	−1.2746	−2.1886	SRC	NM_005417
	−1.2746	−2.9897	CDKN1C	NM_000076
	−1.1892	−3.3173^**	VEGFA	NM_003376
SKBR3	2.2191	1.0792	BAD	NM_004322	Early
	2.639	1.2658	CCNE1	NM_001238
	2	1.0281	IGFBP3	NM_000598
	−2.0705^*	−1.5583	MUC1	NM_001018016

Note: Values are fold changes of ±2 at each time point relative to control with the statistically significant values indicated by asterisks (^*P<0.05 and ^**P<0.01)

Abbreviation: TAM, tamoxifen.

Conversely, T47D was a cell line in which a smaller effect was seen after TAM addition. Only three out of the 84 analyzed genes were regulated (3.6%), and all of them were upregulated only at 24 h (Table 3 and Figure 2B). None of these three genes showed statistically significant differences in expression.

BT474 was the cell line in which TAM showed a higher effect; 11 out of 84 analyzed genes were regulated (13.1%): one gene (9.1%) was upregulated, while the other 10 genes (90.9%) were suppressed (Table 3 and Figure 2C). In addition, the expression of all of these genes was modified only at later time point (48 h). Among these genes, four (36.3%) showed statistically significant altered expression: three were downregulated (EGFR, JUN, and VEGFA), and one (PGR) was upregulated.

Finally, TAM effects on SKBR3 cells were slight, modulating four out of the 84 analyzed genes (4.8%) at 24 h (early expression): three of them were upregulated and one gene was downregulated (Table 3 and Figure 2D). Only the MUC1 gene showed statistically significant downregulation relative to the control.

Pathway analysis of E2-regulated genes

To further evaluate data at the biological level, pathway analysis was conducted by ORA. Table 4 lists biological pathways overrepresented after E2 addition with pathways in which the expression levels of significantly modulated genes were changed with respect to those that would be expected to change by chance.

Table 4 List of biological pathways overrepresented by up- or downregulated genes in MCF7, T47D, BT474, and SKBR3 cells after E2 treatments

Cell line	Time	Signal pathways	Genes	FE	P (Fisher)
Upregulated
MCF7	24 h	Regulation of cell cycle process	BRCA2, BIRC5, EGF, BRCA1	4.7	3.30E−03
	24 h	DNA replication	BRCA2, EGF, BRCA1	7.1	3.50E−03
	48 h	DNA replication	BRCA2, EGF, BRCA1	8.3	2.10E−03
T47D	24 h	Steroid hormone receptor signaling pathway	PGR, CCNE1, ESR2, BRCA1	5.3	2.10E−03
	24 h	Intracellular receptor-mediated signaling pathway	PGR, CCNE1, ESR2, BRCA1	5.3	2.10E−03
	24 h	Regulation of transcription	PGR, CCNE1, BRCA2, ESR2, MYC, BRCA1	2	2.00E−02
	24 h	Cell cycle process	CCNE1, BRCA2, BIRC5, MYC, BRCA1	2.3	2.90E−02
	48 h	Cell cycle process	MKI67, BRCA2, BIRC5, THBS1, MYC, TGFB1	2.7	3.30E−03
	48 h	Cell cycle	MKI67, BRCA2, BIRC5, THBS1, MYC, TGFB1	2.5	6.60E−03
	48 h	Antiapoptosis	BIRC5, ESR2, THBS1, MYC	3.6	1.00E−02
	48 h	Response to abiotic stimulus	BRCA2, THBS1, MYC, TGFB1	3.2	1.80E−02
	48 h	Regulation of mitotic cell cycle	BRCA2, BIRC5, MYC, TGFB1	3.2	1.80E−02
BT474	24 h	Positive regulation of molecular function	NME1, THBS1, MYC	3.9	2.10E−02
	48 h	Antiapoptosis	IGF1R, THBS1, MYC	4.8	1.10E−02
SKBR3	24 h	Developmental growth	CCNE1, BRCA2, PLAU	6.9	3.50E−03
	24 h	Growth	CCNE1, BRCA2, PLAU	5.7	6.60E−03
	48 h	Response to toxin	BCL2, CDH1	27.77	8.80E−04
Downregulated
MCF7	24 h	Regulation of cell development	IGFBP3, SLIT2	6.2	5.50E−03
	48 h	Positive regulation of apoptosis	CDKN1A, MMP9, RARB, NR3C1, IGFBP3	2.2	3.40E−02
	48 h	Positive regulation of cell death	CDKN1A, MMP9, RARB, NR3C1, IGFBP3	2.2	3.40E−02
	48 h	Positive regulation of programmed cell death	CDKN1A, MMP9, RARB, NR3C1, IGFBP3	2.2	3.40E−02

Notes: Pathway analysis was conducted via ORA. ORA was performed in DAVID Bioinformatics Resources 6.7 of the National Institute of Allergy and Infectious Diseases (National Institutes of Health).²⁴ Pathways listed include those whose differentially expressed genes were significantly altered after 24 and 48 h of E2 treatment.

Abbreviations: E2, 17β-estradiol; FE, fold enrichment; ORA, overrepresentation analysis.

In MCF7 cells, E2 stimulated the expression of genes associated with the cell cycle process and DNA replication (BRCA2, BIRC5, EGF, and BRCA1). The BIRC5 gene encoding survivin is a member of the inhibitor of apoptosis gene family that encodes negative regulatory proteins that prevent apoptotic cell death. Amplification of this gene has been reported in 15–30% of BCs, and it has been shown to predict the distant recurrence.²⁶ Similarly, overexpression of BRCA1 and BRCA2 genes can cause an aberrant response to DNA damage. Thus, upregulation of these genes probably leads to an overall increase in both proliferation and cell survival. Conversely, addition of E2 to MCF7 cells suppressed genes involved in regulation of cell development (IGFBP3 and SLIT2), and apoptosis activation and cell death (CDKN1A, MMP9, RARB, NR3C1, and IGFBP3). The IGFBP3 gene has antiproliferative effects on cancer. It is involved in the repair of DNA damage in BC cells²⁷ and controls cell growth by feedback.²⁸ In addition, deregulation or loss of IGFBP3 gene expression has been associated with other carcinomas,²⁹ whereas the SLIT2 gene possesses a suppressive activity against cancer metastasis.³⁰ Deregulation of genes implicated in apoptosis activation and cell death might be related to the increased cell proliferation observed in these cells and could favor tumor development.

In T47D cells, there were several genes upregulated by E2 with known relevant oncogenic functions. Here, we could observe differences in pathway activation according to treatment time points. Thus, at 24 h, E2 activated genes related to signaling pathways of steroid hormones and regulation of transcription (PGR, CCNE1, BRCA2, ESR2, MYC, and BRCA1), while at 48 h, E2 activated genes related to regulation of the cell cycle process and response to abiotic stimulus (MKI67, BRCA2, BIRC5, THBS1, MYC, and TGFB1) (Table 4).

Genes preferentially upregulated by E2 in BT474 cells encode proteins with diverse roles in multiple functional pathways including positive regulation of molecular functions (NME1, THBS1, and MYC) and antiapoptosis (IGF1R, THBS1, and MYC) (Table 4). NME1 is involved in cell proliferation and differentiation, and its high expression has been associated with sporadic colorectal cancer.³¹ The THBS1 gene encodes a glycoprotein that mediates cell-to-cell and cell-to-matrix interactions, and has been shown to play roles in the promotion of angiogenesis and metastasis of invasive ductal BC.³² The MYC gene encodes a protein that promotes cell cycle progression, apoptosis, and cellular transformation. Mutations, overexpression, rearrangement, and translocation of this gene have been observed in various tumors.³³ Elevated expression levels of the IGF1R gene have been associated with poor survival of patients with invasive BC³⁴ and TAM resistance in xenograft models.³⁵

In SKBR3 cells, E2 caused upregulation of genes associated with growth (CCNE1, BRCA2, and PLAU) and responses to toxins (BCL2 and CDH1). These signaling pathways are involved in positive regulation of cell migration, proliferation, and apoptosis, and all of them correlate with carcinogenesis. Of note, BCL2 is one of the most important antiapoptotic genes and widely related to survival of tumor cells. CCNE1 is an important regulator of cell cycle progression, and its overexpression has been associated with poor overall survival and BC-specific survival.³⁶PLAU encodes a serine protease involved in degradation of the extracellular matrix and possibly tumor cell migration and proliferation.³⁷

Notably, in BT474 and SKBR3 cells, genes downregulated by E2 were not found to be significantly associated with any particular signaling pathway. The number of significantly modulated genes by TAM was insufficient to establish any association with a particular signaling pathway in any of the studied cell lines.

Overall representation of the data (PCA)

For a simultaneously general representation of the gene expression response in treated cell lines, PCA was performed (Figure 3). After 24 and 48 h of E2 treatment, MCF7 and BT474 cells had the most diverse conditions, differing in the direction of the positive and negative parts of PC1 compared with the control (time 0). Interestingly, the E2 response was also different in BT474 and MCF7 cells according to the position in PC2 (positive and negative parts). TAM-treated cells had a lesser difference compared with the control (time 0). The largest response after TAM treatment for 24 and 48 h was found in MCF7 cells. However, significant differences were not observed between cell lines. In general, T47D and SKBR3 cell lines had a lesser response in all experimental conditions (Figure 3).

Figure 3 PCA of MCF7, T47D, BT474, and SKBR3 cells after E2 and TAM treatments. PCA was performed for all evaluated genes (84 genes) compared with time 0. Thus grouped, after 24 and 48 h of E2 treatment, MCF7 and BT474 cells showed most diverse conditions compared with the control (time 0) and also between them, differing in the direction of the PC1 and PC2, respectively, while after TAM treatment, the cells had a lesser difference compared with the control (time 0). The largest response after TAM treatment for 24 and 48 h was found in MCF7 cells. Each cell line is represented by a different color, and each treatment performed is represented by a different geometric figure.

Abbreviations: PCA, principal component analysis; E2, 17β-estradiol; TAM, tamoxifen; PC, principal component.

The pcaGoPromoter package²⁵ was used to identify genes responsible for the differences (loading genes) in each PC (Figure 3). PC1 positive describes the lesser response to TAM in most treated cell lines and E2 response in SKBR3 and T47D cells. The top five genes involved in the positive part were CCNE1, BCL2, BAD, PTGS2, and GLI1. PC1 negative and PC2 positive describe the common E2 response in MCF7 and BT474 cells in which the top five genes involved were PGR, MYC, IGFBP3, THBS1, and EGFR. However, the specific E2 response in BT474 cells could be explained by the positive part of PC2, in which the top five genes involved were MYC, IGF1R, THBS1, EGFR, and MUC1 (FC: ±2 and P<0.05; Table 2). Conversely, the specific E2 response in MCF7 cells was explained by the negative part of PC2, in which the top five genes involved were EGF, CSF1, BRCA2, CDKN1C, and IGFBP3 (FC: ±2 and P<0.05; Table 2).

Discussion

One area of great importance in BC research is the study of gene expression regulated by both estrogenic and antiestrogenic agents. Although many studies have been performed in this area, most of them have only examined the effects of E2 and TAM on MCF7 cells.^14,38–41 In the present study, we analyzed changes in gene expression in not only MCF7 cells but also other ERα+ and ERα− cell lines (T47D, BT474, and SKBR3) for 84 genes related to BC development. This allowed us to determine the expression patterns modulated by E2 and TAM in these cells and perform comparisons between other cell lines and MCF7 cells.

After E2 treatment, several patterns of regulation were apparent: while in MCF7 the number of upregulated and downregulated genes was similar, the majority of them were upregulated in T47D cells, and downregulated in SKBR3 and BT474 cells. After E2 addition to MCF7 cells, we obtained similar results to those reported previously in other studies.³⁹ Interestingly, genes such as EGF, RARB, CCND2, IGFBP3, SLIT2, and HPRT1 were significantly regulated only in this cell line, not only indicating specific signal pathway activation but also revealing that HPRT1 would not be useful as a reference/housekeeping gene at least in MCF7 cells. The BIRC5 gene was upregulated in both MCF7 and T47D cells, which confirmed that hormonal regulation of survivin mediated some of the effects of E2 on programmed cell death in ERα+ BC cells.^42,43 Other genes that were upregulated in MCF7 and T47D cells after E2 treatment were MKI67, BRCA1, BRCA2, and PGR, which suggest increased cell responses to genomic damage.³⁹

Although T47D and MCF7 cells are immunohistochemically similar (ERα+/HER2−), we observed differences in regulation of genes such as c-MYC, TGFB1, and THSB1. Nevertheless, although the c-MYC gene has been reported to be upregulated by E2 in MCF7 cells,⁴⁴ we did not find statistically significant changes in its expression. Conversely, such changes were markedly observed in T47D and BT474 cells.

Similar to MCF7 cells, the BT474 cell line (ERα+/HER2+) showed a larger number of changes in gene expression after E2 treatment, and the number of suppressed genes was greater than that of upregulated genes. In addition to the highest upregulation of c-MYC and PGR genes after E2 treatment, BT474 cells showed specific and important regulation of genes such as THBS1, EGFR, IGF1R, and MUC1. Interestingly, genes associated with the HER2 signaling pathway (IGFBP3, PGR, SLIT2, ERBB2, GRB7, and NR3C1) exhibited the same trend in both MCF7 and BT474 cells, although not significant in all cases, indicating a similar response mediated by ER independently of the HER2 status.

SKBR3 cells (ERα−/HER2+) were clearly different from MCF7 cells, because only the BRCA2 gene showed a similar expression pattern after E2 treatment. Interestingly, activation of antiapoptotic BCL2 was specific for SKBR3 cells, and such regulation has been reported as a direct target of transcriptional regulation by ER in MCF7 cells.⁴⁵ Although we did not find regulation of BCL2 in MCF7 cells, these observations strongly support the notion that ERα is not always required for E2 regulation and provide evidence that E2 activates other signaling pathways not associated with ERα. In SKBR3 cells, genes such as IL6, PLAU, and CDKN1C also showed specific and important regulation in at least one of the E2 experimental conditions.

Results from pathway analysis suggested that the gene expression program activated by E2 reflects more than a simple mitogenic response to this hormone in all cell lines. Accordingly, many of the E2-regulated genes were involved in cellular functions other than cell cycle control, such as cell development changes (SLIT2 and IGFBP3) and cell survival (BIRC5).^14,43,46 Although the immunohistochemical status allows supposing a similar response to E2 in ERα+ BC cell lines (MCF7, T47D, and BT474), PCA showed how ERα+ cells behave differently, not only with respect to ERα− cells (SKBR3) but also between them. Accordingly, MCF7 cells exhibited an important activation of DNA replication through genes such as EGF and BRCA2, whereas in BT474 cells, mainly antiapoptotic mechanisms were activated through MYC, IGF1R, and THSB1 genes.

Even if PCA confirmed that gene expression changes after TAM treatment were not significant to differentiate BC cell lines, a smaller effect was observed and some genes displayed similar expression patterns to those seen after E2 treatment. In three of the tested cell lines, both E2 and TAM modulated the expression of PGR and EGF in MCF7 and BT474 cells, and MUC1 in SKBR3 cells. Our results suggest similarities between signaling pathways activated by E2 and TAM. These similarities could be associated with “non-genomic” actions of E2 in target cells, leading to direct hormonal regulation of multiple cytoplasmic signal transductions, as suggested previously.^47–51

Our results are relevant, because many studies have reported that resistance to TAM may be due to its substantial E2 agonist effect in ERα+/HER2+ cells.^{48,49,51–53} Similarly, our data indicate that TAM could perform this E2 agonist role in not only ERα+/HER2+ (BT474) cells but also both ERα+/HER2− (MCF7) and ERα−/HER2+ (SKBR3) cells, according to the cellular response (gene expression pattern) triggered by this drug. In fact, recently it has been indicated that in BC patients, the balance between agonist and antagonist properties of TAM differs among cell types, supporting our observations.⁵⁴

Also, it is important to consider that the effects of E2 and TAM on these cells could be explained by genetic heterogeneity among BC cells and also differences in ERα amount. In fact, it has been reported that T47D cells have low levels of ERα compared with MCF7 cells.⁵⁵ Additionally, it is necessary to note that other ER forms, such as G protein-coupled estrogen receptor (GPER) and ERβ, might also be involved in the effects observed after treatments. GPER, an estrogen transmembrane receptor, which modulates both rapid non-genomic and genomic transcriptional events of estrogens, has been observed in ERα+ cells (MCF7 and T47D).⁵⁶ In addition, recent emerging evidence demonstrated that GPER could be a potential transducer of E2 in SKBR3 cells.⁵⁷ Thus, the presence of this receptor could also explain the modulation of gene expression induced by E2 and TAM in all cell lines analyzed. In fact, GPER has been involved in drug resistance, which often occurs during cancer treatments.⁵⁸ Otherwise, considering that in human BC cell lines the expression of endogenous ERβ has not been clearly described,⁵⁹ it could be suggested that the modulation of gene expression in ERα+ and ERα− cells could be attributed to both ERα and GPER. However, further functional studies are needed to determine if the observed gene expression changes were modulated by ERα, GPER, or both.

Concurrent with our work, some studies identified genes associated with TAM resistance in ERα+ BC patients. Interestingly, some of these genes were found modified by TAM in our study, which supports our approach. One of these studies found a low expression of CDKN1A after TAM treatment,²⁰ as found by us in MCF7 cells. CDKN1A encodes the protein p21^WAF1/CIP1; this protein has been reported to be absent in a clinical case of TAM-stimulated growth.⁶⁰ p21^WAF1/CIP1 interacts with several cell cycle regulators, but the precise mechanism behind its role in TAM resistance is unclear. Additional studies showed significantly lower expression of CDKN1C gene in BC patients treated with adjuvant therapy (TAM)⁶¹ and a high expression of PGR gene in TAM-treated patients.⁶² Similar results were observed by us in ERα+ cell lines (MCF7 and BT474) after TAM addition. CDKN1C is a tumor suppressor gene that regulates cell proliferation and has recently been found to be downregulated in metastatic tumors. PGR gene encodes a protein that mediates the physiological effects of progesterone, which plays a central role in reproductive events associated with the establishment and maintenance of pregnancy. Interestingly, this gene has recently been considered as one of the most promising prognostic biomarkers in TAM-treated patients.⁶²

Taking all together, our results could provide important markers to be validated in BC patient samples, and subsequently used for predicting the outcome in ERα+ and ERα− tumors after TAM or hormonal therapy. Exploring novel markers for TAM therapy will be helpful for the treatment of BC.

Conclusion

Although the MCF7 cell line is the most frequently used model to determine gene expression profiles in response to treatments with E2 and TAM, the changes in gene expression patterns identified in ERα+ and ERα− cell lines could reflect distinctive properties of these cells, which could be exploited not only to identify the response of various cell types to E2 but also to examine whether these genes might serve as markers of TAM sensitivity and/or participate in the development of TAM resistance. Indeed, we have shown that the expression patterns of some genes can be used to discriminate between ERα+ cells, and between ERα+ and ERα− cells.

Author contributions

Milena Rondón-Lagos and Nelson Rangel conceived and designed the study, and wrote the manuscript. Victoria E Villegas contributed to writing of the manuscript and made critical revisions. All authors contributed toward data analysis, drafting and critically revising the paper and agree to be accountable for all aspects of the work.

Disclosure

The authors report no conflicts of interest in this work.

Acknowledgments

Nelson Rangel is supported by a Colciencias Grant (call 617, Colombia). The authors acknowledge the scientific support provided by Professor Anna Sapino. César Payán-Gómez is acknowledged for bioinformatic analyses. The authors thank Mitchell Arico from Edanz Group (www.edanzediting.com/ac) for editing the draft of this manuscript.