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,[29]^1 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.[30]^2 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.[31]^3 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,[32]^4 through transcription factor cross talk,[33]^5 and non-genomic mechanisms whereby E2 binds to ERs localized in the cell membrane and activates signal transduction pathways.[34]^6 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.[35]^7^,[36]^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.[37]^9^–[38]^11 However, responses to TAM have frequently limited durations due to the development of resistance.[39]^12^,[40]^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.[41]^14^,[42]^15 Considering these aspects, several studies have evaluated gene expression profiles, which focused on identifying changes in gene expression associated with BC prognoses,[43]^16^–[44]^19 drug resistance,[45]^20 or tumor aggressiveness.[46]^21 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[2]. 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.[47]^22 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^2 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^2 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^2 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 ([48]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^2 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 [49]Open in a new tab 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.[50]^23 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)[51]^24 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.[52]^23 To identify the main genes responsible for the different PC segments, we used the pcaGoPromoter package[53]^25 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[54]^25 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 ([55]Table 2 and [56]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) ([57]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 ([58]Table 2). 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 [59]NM_001168 Early 2.1435^* 1.7411 NME1 [60]NM_000269 −2.4623 −1.2311 SERPINE1 [61]NM_000602 2.4623^** 2.639^** BRCA1 [62]NM_007294 Early and late 3.249^** 2.3784^** BRCA2 [63]NM_000059 −2^** −2.0705^** CCND2 [64]NM_001759 3.4822^* 5.0982^** EGF [65]NM_001963 2.8284^** 2.4623^** HPRT1 [66]NM_000194 −3.0314 −2.0705 TWIST1 [67]NM_000474 −2.4623 −2.4623 ERBB2 [68]NM_004448 −2.8945 −3.1383 GRB7 [69]NM_005310 −2.639^** −5.0982^** IGFBP3 [70]NM_000598 3.249^** 2.0705^** MKI67 [71]NM_002417 6.0629^** 8.8766^** PGR [72]NM_000926 −2.1435^** −2.2191^** RARB [73]NM_000965 −2.639^** −3.7321^** SLIT2 [74]NM_004787 1.7613 2.5491^** TFF3 [75]NM_003226 Late −1.7211 −2.1435^** NR3C1 [76]NM_000176 −1.7211 −2.2974^** CDKN1A [77]NM_000389 T47D 2.9794^** 1.9656 BRCA1 [78]NM_007294 Early 2.7678^* 1.9571 CCNE1 [79]NM_001238 2.1067^** 1.7715 PGR [80]NM_000926 2.9409 1.8742 TP73 [81]NM_005427 2.6735 2.3274 ABCG2 [82]NM_004827 Early and late 2.0885^* 2.4837^* BIRC5 [83]NM_001168 2.3784^** 5.278^** MYC [84]NM_002467 2.2482^** 2.4095^** BRCA2 [85]NM_000059 2.3784 2.3784 SFN [86]NM_006142 2.2875 2.2095 SRC [87]NM_005417 2.7203 3.716^** ESR2 [88]NM_001437 1.858 2.7203^** MKI67 [89]NM_002417 Late 1.5223 2.0795 NOTCH1 [90]NM_017617 1.9487 2.0885 RASSF1 [91]NM_007182 1.3779 3.3782^* TGFB1 [92]NM_000660 1.3899 2.2579^** THBS1 [93]NM_003246 1.5422 2.2776 CTSD [94]NM_001909 BT474 −2.5669^* −1.3013 SNAI2 [95]NM_003068 Early 2.2815^** 1.8277 NME1 [96]NM_000269 3.1167^** 3.6553^** THBS1 [97]NM_003246 Early and late −2.5669 −2.7895 GATA3 [98]NM_002051 3.1167^** 2.8679^** IGF1R [99]NM_000875 −2.4794^** −3.9449^** EGFR [100]NM_005228 −4.9588 −2.8879 IGFBP3 [101]NM_000598 25.8125^** 15.1369^** MYC [102]NM_002467 9.1261^** 9.3179^** PGR [103]NM_000926 −2.3134^* −2.6945^* MUC1 [104]NM_001018016 −2.2346 −2.042 SLIT2 [105]NM_004787 −1.7532 −2.514 ERBB2 [106]NM_004448 Late −1.3287 −2.114 ESR1 [107]NM_000125 −1.4743 −2.2658 FOXA1 [108]NM_004496 −1.5263 −2.114^* GRB7 [109]NM_005310 −1.0792 −2.042^* GSTP1 [110]NM_000852 −1.815 −2.042^** NR3C1 [111]NM_000176 1.3566 4.0558^** ABCG2 [112]NM_004827 SKBR3 2.0849 1.1096 BAD [113]NM_004322 Early 2.0849^* 1.1096 BRCA2 [114]NM_000059 3.7581^** 1.366 CCNE1 [115]NM_001238 3.2716 1.366 IGFBP3 [116]NM_000598 5.1337^** 1.4142 IL6 [117]NM_000600 −2.4453^** −1.2746 MUC1 [118]NM_001018016 2.8481^* 1.9319 PLAU [119]NM_002658 −2.1287 −1.0718 PYCARD [120]NM_013258 2.395 1.1096 SERPINE1 [121]NM_000602 −2.1287^* −2.2974^** CDKN1C [122]NM_000076 Early and late 4.0278 2.1435^** BCL2 [123]NM_000633 1.3566 2.2191^** CDH1 [124]NM_004360 Late [125]Open in a new tab 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. Figure 1. [126]Figure 1 [127]Open in a new tab 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. 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 ([128]Table 2 and [129]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 ([130]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 ([131]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 ([132]Table 2 and [133]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 ([134]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 ([135]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 ([136]Table 2 and [137]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 ([138]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 ([139]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 ([140]Table 3 and [141]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). 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 [142]NM_000076 Early −2.0705 −1.6434 CSF1 [143]NM_000757 3.1383^* 1.977 PGR [144]NM_000926 4.2871^** 3.6893^** EGF [145]NM_001963 Early and late −1.366 −2.0232 SERPINE1 [146]NM_000602 Late T47D 2.3681 1.9235 GATA3 [147]NM_002051 Early 2.4516 1.3601 RARB [148]NM_000965 2.2677 1.2152 TWIST1 [149]NM_000474 BT474 −1.2746 −2.4284^** EGFR [150]NM_005228 Late −1.5157 −2.6027 ERBB2 [151]NM_004448 −1.7411 −2.3457 GATA3 [152]NM_002051 −1.1892 −2.114 GRB7 [153]NM_005310 −1.1096 −2.114^** JUN [154]NM_002228 −1.2311 −2.3457 NOTCH1 [155]NM_017617 1.8025 2.6759^* PGR [156]NM_000926 −1.4641 −2.6945 RASSF1 [157]NM_007182 −1.2746 −2.1886 SRC [158]NM_005417 −1.2746 −2.9897 CDKN1C [159]NM_000076 −1.1892 −3.3173^** VEGFA [160]NM_003376 SKBR3 2.2191 1.0792 BAD [161]NM_004322 Early 2.639 1.2658 CCNE1 [162]NM_001238 2 1.0281 IGFBP3 [163]NM_000598 −2.0705^* −1.5583 MUC1 [164]NM_001018016 [165]Open in a new tab 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. Figure 2. [166]Figure 2 [167]Open in a new tab 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. 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 ([168]Table 3 and [169]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 ([170]Table 3 and [171]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 ([172]Table 3 and [173]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. [174]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 [175]Open in a new tab 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).[176]^24 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.[177]^26 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[178]^27 and controls cell growth by feedback.[179]^28 In addition, deregulation or loss of IGFBP3 gene expression has been associated with other carcinomas,[180]^29 whereas the SLIT2 gene possesses a suppressive activity against cancer metastasis.[181]^30 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) ([182]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) ([183]Table 4). NME1 is involved in cell proliferation and differentiation, and its high expression has been associated with sporadic colorectal cancer.[184]^31 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.[185]^32 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.[186]^33 Elevated expression levels of the IGF1R gene have been associated with poor survival of patients with invasive BC[187]^34 and TAM resistance in xenograft models.[188]^35 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.[189]^36 PLAU encodes a serine protease involved in degradation of the extracellular matrix and possibly tumor cell migration and proliferation.[190]^37 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 ([191]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 ([192]Figure 3). Figure 3. Figure 3 [193]Open in a new tab 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[194]^25 was used to identify genes responsible for the differences (loading genes) in each PC ([195]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; [196]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; [197]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.[198]^14^,[199]^38^–[200]^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.[201]^39 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.[202]^42^,[203]^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.[204]^39 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,[205]^44 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.[206]^45 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).[207]^14^,[208]^43^,[209]^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.[210]^47^–[211]^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.[212]^48^,[213]^49^,[214]^51^–[215]^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.[216]^54 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.[217]^55 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).[218]^56 In addition, recent emerging evidence demonstrated that GPER could be a potential transducer of E2 in SKBR3 cells.[219]^57 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.[220]^58 Otherwise, considering that in human BC cell lines the expression of endogenous ERβ has not been clearly described,[221]^59 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,[222]^20 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.[223]^60 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)[224]^61 and a high expression of PGR gene in TAM-treated patients.[225]^62 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.[226]^62 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. Acknowledgments