An ERRα-ZEB1 transcriptional signature predicts survival in triple-negative breast cancers

doi:10.21203/rs.3.rs-4869822/v1

Download PDF

Research Article

An ERRα-ZEB1 transcriptional signature predicts survival in triple-negative breast cancers

https://doi.org/10.21203/rs.3.rs-4869822/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background.

Transcription factors (TFs) act together with co-regulators to modulate the expression of their target genes, which eventually dictates their pathophysiological effects. Depending on the co-regulator, TFs can exert different activities. The Estrogen Related Receptor α (ERRα) acts as a transcription factor that regulates several pathophysiological phenomena. In particular, interactions with PGC-1 co-activators are responsible for the metabolic activities of ERRα. In breast cancers, ERRα exerts several tumor-promoting, metabolism-unrelated activities that do not depend on PGC1, questioning the identity of the co-activators involved in these cancer-related effects.

Methods.

We used bio-computing methods to identify potential co-factors that could be responsible for the activities of ERRα in cancer progression. Experimental validations were conducted in different breast cancer cell lines, using determination of mRNA expression, ChIP-qPCR and proximity ligation assays.

Results.

ZEB1 is proposed as a major ERRα co-factor that could be responsible for the expression of direct ERRα targets in triple-negative breast cancers (TNBC). We establish that ERRα and ZEB1 interact together and are bound to the promoters of their target genes that they transcriptionally regulate. Our further analyses show that the ERRα-ZEB1 downstream signature can predict the survival of the TNBC patients.

Conclusions.

The ERRα-ZEB1 complex is a major actor in breast cancer progression and expression of its downstream transcriptional targets can predict the overall survival of triple-negative breast cancer patients.

transcription

triple-negative breast cancer

co-regulators

ERRα

ZEB1

Transcription factors (TFs) can be defined as proteins that bind to discrete response elements on the promoter and/or enhancer of their target genes and regulate their expression. However, to achieve the latter function, TFs require the binding of additional co-factors. Directly or indirectly, these factors will induce modifications of the chromatin on the vicinity of the target transcriptional start site, allowing the recruitment of active or repressive transcription modulating complexes. Literature shows an increased interest in studying the functional interactions between transcription factors and their co-factors. For instance, members of the TEAD family can recruit the YAP/TAZ co-activators to enhance the expression of their targets [1]. In contrast, binding of the VGLL4 co-repressor to TEAD will reduce the expression of these targets [2]. Though the Glucocorticoid Receptor (GR) acts as a transcriptional activator on its corresponding binding sites, it will antagonize the activation exerted by AP1 or NFκB on genes involved in inflammation (reviewed in [3]). Co-factors thus appear as instrumental in determining the activities of individual TFs.

ERRα (encoded by ESRRA) is a member of the nuclear receptor superfamily that acts as a transcription factor involved in various pathophysiologic effects. For instance, this protein is highly active in the regulation of cellular metabolism that it regulates by binding to the promoters of its metabolic targets such as MCAD, IDH3A and ATP5B and recruitment of members of the PGC-1 family of co-activators (reviewed in [4, 5]). These functions are operated in tissues displaying a high energy demand such as liver, muscle or adipose. Remarkably, in breast cancer cells, ERRα does not regulate these metabolic targets although it is bound to their promoters [6, 7]. In contrast, ERRα displays an array of cancer-related functions in these cells through the regulation of distinct transcriptional targets, such as WNT11, TNFAIP1 or MMP1 (reviewed in [8]). The receptor indeed modulates proliferation, resistance to hypoxia, angiogenesis, cell adhesion and migration, as well as extracellular matrix invasion [6, 9–13]. It is thus possible that co-factors differing from PGC-1s and that remain to be identified contribute to these activities. These data also suggest that ERRα could be used as a pharmacological target to be inactivated in various types of cancers. However, inhibiting the receptor could also impact its other activities, and could therefore yield unwanted and detrimental side effects. As an example, mice depleted from ERRα do not resist to cold exposure [14]. Identifying ERRα co-factors that are selectively involved in its cancer functions could suggest new targets to be examined against cancers. Several factors have already been identified as interacting with ERRα and modulating its cancer-promoting functions. This is the case for NCOA3 (AIB1, ref. [15]), NCOA2 (SRC-2, ref. [16]), KDM1A (LSD1, ref. [12]), SETD7 [7] or NRIP1 (RIP140, ref. [17]). These factors have been identified on a gene candidate basis and this illustrates the lack of a general unbiased method to identify ERRα coregulators outside metabolic cells.

Here we used direct ERRα transcriptional targets that we have identified by RNA-Seq and ChIP-Seq in triple-negative breast cancer (TNBC) cells. This gene list was submitted to bio-computing approaches to suggest transcription factors that, together with ERRα, could be involved in their high expression in breast cancer samples. These methods suggested ZEB1 as the main factor responsible for these activities. This hypothesis was experimentally verified. We focus on eight target genes that we show to be transcriptionally regulated by ZEB1 and ERRα in TNBC cells. These two factors interact together and are recruited to ERRα response elements on the promoter of their target genes. Importantly, high expression of their transcriptional signature confer a poor survival prognosis to TNBC.

Cells. MCF7 (ER+/PR+/Her2-), BT474 (ER+/PR+/Her2+), SKBr3 (ER-/PR-/Her2+) and MDA-MB231 (TNBC; mesenchymal-like subtype) cells were cultured in DMEM supplemented with 10% fetal calf serum (FCS), 10 U/ml penicillin and 10 µg/ml streptomycin. BT549 (TNBC; mesenchymal subtype), HCC38 (TNBC; basal-like 1 subtype) and MDA-MB468 (TNBC; basal-like 1 subtype) cells were cultured in RPMI supplemented with 10% FCS, 10 U/ml penicillin and 10 µg/ml streptomycin. For siRNA transfection, 3.10^− 5 cells per ml were seeded in 6-well plate and 25 pmol/ml of total siRNA were transfected with INTERFERin (Polyplus Transfection; Illkirch-Graffenstaden, France) according to the manufacturer’s recommendations. siRNAs were from Invitrogen (Strasbourg, France) and sequences are shown on Additional File 2, Table S1.

Expression analyses. Total RNAs were extracted by the guanidinium thiocyanate/phenol/chloroform method. 1 µg of RNA was converted to first strand cDNA using the RevertAid kit (ThermoScientific, Strasbourg, France). Real time PCRs were performed in 96 well plates using the IQ SYBR Green Supermix (Bio-Rad; Marnes-la-Coquette, France). Data were quantified by the ΔΔ-Ct method and normalized to 36b4 expression. Sequences of the primers used in this study are shown on Supplementary Table S1. For western blot assays, cells were lysed in RIPA buffer supplemented with protease inhibitor cocktail (Sigma-Aldrich, Darmstadt, Germany). Proteins (25 µg) were resolved on 7 to 10% SDS-PAGE, blotted onto PVDF membrane (Bio-Rad) and probed with specific antibodies after saturation. The antibodies used in this study were: GAPDH (G9545, Sigma-Aldrich,1/1000), ERRα (GTX108166, Genetex, 1/2500), ZEB1 (HPA027524, Sigma-Aldrich, 1/1000), TFAP2C (sc-12762, 1/1000, Santa Cruz, Heidelberg, Germany). For proximity ligation assays cells cultured on coverslips were fixed with 2% paraformaldehyde, washed with PBS, and analyzed with Duolink PLA kit (Sigma-Aldrich) according to the recommendations from the manufacturer. Samples were DAPI-counterstained. Images were acquired using a Zeiss AxioImager microscope. Analyses were performed as previously described [54].

ERRα directly activated genes. For RNA-sequencing, RNAs were extracted from MDA-MB231 cells cultured in suspension after siRNA treatment, using Qiagen extraction kit (Courtaboeuf, France). RNA-seq was performed on triplicate samples and libraries built using the mRNA-Seq Library Prep kit of Lexogen following the instructions for 5500 SOLiD and including a conversion step to SOLiD 5500W. Sequencing was performed with SOLiD 5500W System (Life Technologies, Strasbourg, France). Sequences were aligned on the human genome (hg38 version) by STAR (v2.7.3a) tool. Read counts were determined using HTSeq_count (Python v3.8.3). Differentially expressed genes were identified with R package DESeq2 (v1.40.2) in R software (v4.3.1), using an adjusted p-value < 0.05 and fold change threshold of 1.5 and 0.75 for over- and under-expressed genes respectively. Genes showing significantly modified expression with both siERRα were considered as ERRα-modulated genes and referred to as differentially expressed genes (DEGs). This list was intersected with the 4846 genes displaying ERRα peak in ChIP-Seq [7] which yielded 190 ERRα directly activated DEGs (Additional File 2, Table S2). GO analysis was performed on these 190 DEGs by R package goseq (v1.52.0). Heatmaps were generated using R package pheatmap (v1.0.12). Bar plots were generated using R package ggplot2 (v3.4.4). Results of the RNA-Seq are available at GEO website under accession number GSE198923.

Transcriptional Regulators (TRs) collection. 2175 human TRs were collected from public databases as in Cerutti et al. [7]. TRs with mean count < 10 in TCGA breast cancers were removed as well as those with mean expression < 0.1 or standard deviation < 0.1 in upper-quartile normalized data. Log2-transformation was performed on upper-quartile normalized expression and this list was used in PCA analysis by R package FactoMineR (v2.9). TRs with PCA-component 1 or 2 > 0.6 were retained. 38 additional TRs from Cerutti et al., [7] were added producing a final list of 493 pre-selected breast cancer related TRs (Additional File 2, Table S3).

Public expression datasets. TCGA breast cancer RNA-Seq data and clinical records from 1091 primary breast tumors were downloaded from the Genomic Data Commons (GDC) data portal (https://portal.gdc.cancer.gov/). RNA-Seq data and clinical features from 360 Fudan University Shanghai Cancer Center (FUSCC) TNBC were downloaded from ref. [19]. Breast cancer expression data, gene annotation and clinical features were downloaded for 180 primary breast tumors from Gene Expression Omnibus (GEO) database (accession number GSE18229). 2509 breast cancer mRNA expression z-scores (log microarray) and clinical features in METABRIC were downloaded from https://www.cbioportal.org. RNA-Seq data for 56 breast cancer cell lines from the Cancer Cell Line Encyclopedia (CCLE) were downloaded from DepMap Public 22Q4 (https://depmap.org/portal/download/). scRNA-Seq data from breast cancer cells were collected from ref. [20]. A total 35276 individual cells from 32 human breast cell lines (31 cancer cell lines and 1 non-cancer cell line) were sequenced and recorded by cancer subtypes. Cell lines were subtyped according to ref. [21–23]. Tumors were subtyped according to ref. [23]. Tumors were also sorted according to predicted menopausal status based on WHO recommendations (https://www.who.int/news-room/fact-sheets/detail/menopause) and [55].

Sparse PLS models. An adaptive sparse Partial Least Square (sPLS) regression method in R package plsgenomics (v1.5-2) was used as described in ref. [7]. Briefly, the value of latent components and lambda parameter were firstly tuned and determined by 10-fold cross-validation. After expression modeling based on these optimized parameters, R-squared determination coefficient (R²) was used to evaluate each gene model’s quality. Gene expression models described by a set of TRs in high-quality (R² > 0.6) were selected and kept for TR prediction. For each DEG, computation was run 10 replicates. Contribution of each TR to the DEGs expression (mean-prop) was counted as the number of genes including this TR in their high-quality computation models. TRs with the highest mean-prop values were classed as best predicted ERRα co-factor. Modeling was performed three times (i.e. 30 replicates altogether) to evaluate the robustness of the sPLS method. Predicted results were displayed in circular bar plot by using R package ggplot2 (v3.4.4).

Gene expression comparison between breast cancer subtypes. Gene expression comparison in the breast cancer single cell atlas was obtained using http://bcatlas.tigem.it [20]. Gene expression comparison in human breast cancer cell lines and tumors were based on normalized expression data in log-transformed in R. For each gene, all samples’ expression were normalized to [-1,1] by:

$$\:{X}_{norm}=\:\frac{2\text{*}(X-{X}_{min})}{{X}_{max}-{X}_{min}}\:\--\:1$$

and then were compared between molecular subtypes. Expression heatmaps were generated using package pheatmap (v1.0.12).

Gene expression and EMT score. EMT scores were quantified in each breast cancer tumor sample [56]. Briefly, 144 epithelia (Epi) and 170 mesenchymal (Mes) signatures were used to calculate their cumulative distribution function (CDF). EMT score was obtained based on the maximum distance between CDFs of Epi- and Mes signatures using two-sample Kolmogorov-Smirnov test. Pearson correlation test was used to correlate the EMT score to log-transformed gene expression. Regression plots were generated using R package ggplot2 package (v3.4.4). Correlation heatmaps were generated using R package ellipse (v0.5.0).

Survival analyses. Survival analysis in the TCGA breast cancer RNA-Seq dataset was performed using the R package survival (v3.5-7). Each tumor sample was labeled as high or low for each gene as in ref. [23]. High sample is characterized by gene expression (upper-quartile normalized and then log-transformed) higher than the 60th percentile. Low sample displays gene expression lower than the 40th percentile. For gene set combinations, the high/low group comprises samples containing at least proper number (60% * the size of the combination set, for example, when we explore the survival outcome prediction probability of 8 DEGs together, the proper number is 60% * 8 = 5) of high/low labels except in two gene analysis in which proper number was set to be 2. All Kaplan-Meier curves were generated by statistical comparison between two groups based on log-rank test using R package survminer (v0.4.9). Survival analysis of GSE96058 breast cancer RNA-seq dataset was performed by KM Plotter database (https://kmplot.com/analysis/) [57]. For each gene, samples were divided into three terciles based on gene expression. Highest and lowest expression groups were used for further analyses. For the gene set combinations, mean expression of all genes was calculated and used to classify samples to high or low group using the same method in single gene above. Survival analysis in each breast cancer subtype was also conducted by restricting patients to different subtypes (Basal, LuminalA, LuminalB and Her2) in PAM50 subtype selection module. For the breast cancer chip datasets, the following microarray probes were used here: 221421_s_at (ADAMTS12), 204677_at (CDH5), 225442_at (DDR2), 214447_at (ETS1), 225320_at (MCU), 208983_s_at (PECAM1), 206470_at (PLXNC1) and 200665_s_at (SPARC).

ChIP. 10⁷ cells were cross-linked with 1% formaldehyde and quenched with 1 M glycine. After centrifugation, chromatin was extracted from cell pellets and immunoprecipitated using the iDEAL ChIP kit (Diagenode, Liège Belgium) according to the manufacturer’s instructions. Quantitative PCRs were performed using 2 µl of DNA in duplicate and enrichment was calculated related to input.

We investigated the mechanisms through which ERRα positively regulates the expression of its target genes in breast cancers (workflow presented on Additional File 1, Figure S1). We first performed a transcriptomic analysis of this factor in a cellular model (Fig. 1a). To this end, MDA-MB231 cells were grown in culture medium droplets after transfection with two different siRNAs directed against ERRα and RNA sequencing was performed. We observed that the expression of 917 genes (referred to as DEGs; Differentially Expressed Genes) were strongly deregulated by both siRNAs (Additional File 2, Table S2). These included 629 genes that displayed a reduced expression in the absence of ERRα (i.e. stimulated by the receptor). Comparison with our previously published ChIP-Seq data [7] revealed that 190 of these genes presented an enrichment peak for ERRα with a consensus ERRE (ERRα response element) at peak summit. We thus focused on these directly activated DEGs (listed on Additional File 2, Table S2). Gene Ontology (GO) analysis of this gene list showed a massive enrichment in terms related to cell migration (Fig. 1b), consistent with the demonstrated role of ERRα in this phenomenon [6, 10, 18]. This GO term comprised 45 DEGs on which further analyses were conducted. We next established a list of transcriptional regulators (TRs) expressed in breast cancer (Additional File 1, Figure S2a). To this end, a comprehensive set of 2175 TRs was first collected from different databases [7] and submitted to expression-based filtering and PCA analysis of RNA-Seq data from TCGA breast cancers. Altogether, this produced a list of 493 breast cancer relevant TRs (Additional File 2, Table S3) that were used for model computation.

To identify potential ERRα co-activators, we used sparse Partial Least Square (sPLS) regression method for the 45 individual DEGs and the 493 TRs selected above. Modeling was first performed in triplicate on the breast cancer data from TCGA (Additional File 1, Figure S2b), using all 1091 breast primary tumors as well as each subtype in which the data were distributed (luminal-like, Her2+, TNBC or others) (Fig. 2a). Potential co-activators were then ranked and the five top ones per category are shown on Fig. 2a. A certain level of variability was observed within the three replicates and from one subtype to the other. However, for all computed models, ZEB1 was constantly ranked as the top 1 potential co-activator. In the TCGA dataset, the TNBC category only comprises 115 samples (Additional File 1; Figure S2b). To increase the number of TNBC in our study, we used the Fudan University Shanghai Cancer Center (FUSCC) tumor dataset that comprises 360 TNBCs, distributed into four distinct categories (Additional File 1, Figure S2c; [19]). Using the same modelling approach, similar results were obtained that constantly indicated ZEB1 as the top potential co-activator in all TNBC from FUSCC as well as in TNBC subtypes (Fig. 2a). Moreover, we found that 8 (ADAMTS12, CDH5, DDR2, ETS1, MCU, PECAM1, PLXNC1 and SPARC) out of 45 DEGs are particularly sensitive to ZEB1-induced regulation in TCGA tumors (Fig. 2b). This is true for all categories of tumors, including TNBC. We also analyzed the FUSCC TNBC tumors and found that the 8 DEGs are also strongly regulated by ZEB1 in these tumors as well as in their different sub-categories. Modeling on breast tumor datasets that were sorted according to other parameters did not allow to discriminate the ethnic origin of the patients, or their likely pre- or post-menopausal status (Additional File 1, Figure S3a and S3b). Altogether, our data suggest 8 DEGs as reliable ERRα-ZEB1 targets in breast cancer tumors.

To explore the transcription regulation ability of ERRα and ZEB1 together in breast cancer, we then investigated the expression level of ESRRA (which encodes ERRα), ZEB1 and their 8 common targets evidenced above in human breast cancer single cells, cell lines and tumors. Human breast cancer single cell RNA-Seq data (scRNA-Seq) of 35,276 cells from 32 human breast cell lines were collected, processed and analyzed by Gambardella et al. [20]. Six cell clusters (luminal A, luminal B, Her2-enriched, triple-negative type A, triple-negative type B and basal-like) annotated with breast cancer subtypes were identified and showed in an atlas (Fig. 3a). In this atlas, the expression of ESRRA was scattered across all subtypes while ZEB1 expression was enriched in triple-negative subtypes. Interestingly, we observed that the expression of half of 8 DEGs (ADAMTS12, DDR2, ETS1, SPARC) also displayed an enrichment in TNBC subtypes similar to that of ZEB1 (Fig. 3a and Additional File 1 Figure S3c). Next, we investigated RNA expression data of 56 breast cancer cell lines with molecular subtype classification from Cancer Cell Line Encyclopedia (CCLE) [21–23]. We observed that ZEB1 and most of 8 DEGs (ADAMTS12, DDR2, ETS1, MCU, SPARC) have higher normalized expression in the aggressive claudin-low subtype, a sub-fraction of TNBC, while ESRRA displayed no obvious expression preference (Fig. 3b). Four ZEB1-repressed targets (CDH1, LLGL2, CLDN7, CLDN3) were used as control and showed opposite expression pattern to ZEB1 which confirmed the specificity of ZEB1 expression and activity. In addition, we checked their mRNA expression in human breast tumors with identical molecular subtype annotation in cell lines [23]. As expected, claudin-low subtype had an elevated expression of ZEB1 and most of 8 DEGs (CDH5, DDR2, ETS1, PECAM1, PLXNC1, SPARC) (Fig. 3c). Again, the four ZEB1-repressed targets displayed lower expression. In summary, expression profiles demonstrated the higher expression of ZEB1 and most of 8 DEGs in the aggressive claudin-low subtype which strongly suggested ZEB1 could be a co-activator of ERRα.

To experimentally evaluate the effects of ERRα and ZEB1 on the expression of their potential target genes, MDA-MB231 cells were transfected with siRNAs targeting each factor (Fig. 4a). The mRNA expression of the eight DEGs defined above was then studied by RT-qPCR. Depletion of ERRα (Fig. 4b) or ZEB1 (Fig. 4c) reduced the expression of all these genes with the exception of DDR2. Strikingly, simultaneous depletion of both ERRα and ZEB1 did not further decrease target gene expression relative to the inactivation of a single factor (Fig. 4d). This suggests that ERRα and ZEB1 act on the DEGs through the same molecular pathway.

We then analyze the extent of ERRα-ZEB1 cooperation in other breast cancer cell lines. A preliminary RT-qPCR screen showed that, in contrast to ERRα, ZEB1 was weakly expressed (if at all) in three non-TNBC cell lines (MCF7, SKBr3 and BT474) (Ct > 30; Fig. 4e). Interestingly, ERRα depletion in these cells did not globally impact the expression of the DEGs (Additional File 1, Figure S4a). This is not due to a general impairment of ERRα activities in these cell lines, since the receptor is active on other of its transcriptional targets such as LMNB1 or NUDT19. In contrast, ZEB1 and ERRα were strongly expressed in all tested TNBC cell lines (MDA-MB231, HCC38, MDA-MB468 and BT-549; Fig. 4e). Consistently, siRNA-mediated depletion of ERRα or ZEB1 led to reduced expression of all DEGs expressed in these cell lines (Fig. 4b-c and Additional File 1, Figure S4b). Altogether, these data suggest that ERRα-ZEB1 regulate the expression of the eight DEGs specifically in TNBC cell lines.

The extent of these cooperative effects was next evaluated in TNBC cells. To this end, we first analyzed the expression of four ZEB1 targets [23]. Depending on the cell type, CCN2 (CTGF), CCN1 (CYR61) and FOSL1, but not YAP1, can be transactivated by both ZEB1 and ERRα (Additional File 1, Figure S5a). In contrast, ZEB1 does not regulate the expression of GOT2, NOP2, SINHCAF or TMEM120B (Additional File 1, Figure S5b), which are other demonstrated ERRα targets [7]. This indicates that ERRα-ZEB1 cooperation is restricted to a number of genes that do not include all ERRα targets but may extend to several other ZEB1 ones. The effect of ZEB1 is specific of this co-regulator. Indeed, siRNA-mediated depletion of PGC-1α (encoded by PPARGC1A), another ERRα coactivator [24], does not impact on the expression of the 8 DEGs (Additional File 1, Figure S5c).

We next searched to determine the mechanisms through which ERRα and ZEB1 regulate the expression of the DEGs. Analysis of our previously performed ChIP-Seq experiment [7] had shown the binding of ERRα on discrete sites (ERREs; ERRα response elements; 5’-TCAAGGTCA-3’) of the DEGs. This binding was confirmed by independent ChIP-qPCR experiments (Fig. 5a). Strikingly, ChIP-qPCR also demonstrated ZEB1 binding to the same ERRα-binding regions (Fig. 5b), although no conventional direct ZEB1 binding sites (5’-CAGGTG/A-3’) were detected. Binding of ERRα, but not of ZEB1, was also shown on GOT2, NOP2, SINHCAF and TMEM120B promoters (Additional File 1, Figure S5d) in consistence with their response to ERRα but not to ZEB1. This suggests that binding of ZEB1 is limited to a number of ERRα binding sites.

Enrichment of ERRα and ZEB1 in the same genomic regions suggests that the DNA binding of these factors may depend on each other. To examine this possibility, we performed ChIP-qPCR in cells after inactivation of one of the factors. We noted that, in the absence of ERRα, ZEB1 binding to the majority of all DEGs was considerably reduced (Fig. 5c). In contrast, ERRα binding on its corresponding sites was unaffected by ZEB1 depletion (Fig. 5d). Altogether, this shows that ERRα directly binds to its cognate response elements, whereas ZEB1 indirectly contacts these sequences through ERRα. This also suggests that ERRα and ZEB1 interact with each other, a hypothesis that we tested by proximity ligation assay (PLA). Our data indeed document these interactions that took place in the cell nuclei and were impaired upon depletion of ERRα or ZEB1 (Fig. 5e).

ZEB1 is an instrumental factor driving the epithelium-to-mesenchyme transition (EMT) [25–27]. This suggests that its target genes, in particular the DEGs that it shares with ERRα, may be involved in the EMT. To examine this hypothesis, we evaluated the EMT score of each of 1091 breast cancers reported in TCGA and linked this score to the level of expression of ERRα, ZEB1 and each of the DEGs (Fig. 6a). The EMT scores were correlated with the expression level of ZEB1, as expected, but not with that of ERRα. Strikingly, all DEGs except MCU also displayed a positive correlation with the EMT score. We next examined these correlations in three other breast cancer datasets (Fig. 6b). In all these datasets, expression of ZEB1 and most of the examined ZEB1-only (i.e. not regulated by ERRα) targets was correlated to the EMT score. In contrast, expression of ERRα and of ERRα-only (i.e. not regulated by ZEB1) targets was not linked to the EMT score, suggesting that ERRα is not involved in EMT. However, we observed that expression of all ERRα-ZEB1 targets (again except MCU) was correlated to the EMT score, suggesting that the ZEB1-dependent ERRα activity is involved in EMT.

EMT favors the establishment of metastasis and thus decreases survival [28]. The possible involvement of ERRα targets in EMT suggests that these genes may impact patient survival. To test this hypothesis, we examined the overall survival (OS) of breast tumor patients as a function of high or low gene expression. Expression of ERRα or ZEB1, examined individually or together, did not impact on OS (Additional File 1, Figure S6a). In combination, expression of the 8 DEGs did not influence OS. This phenomenon was observed when considering all tumors and also luminal A, luminal B and Her2 + tumors. Strikingly in TNBC, high combined expression of the 8 DEGs reduced OS (Fig. 6c). Again, expression of individual or combined ERRα or ZEB1 did not impact on OS. We next examined the effect of high combined expression of the 8 DEGs on OS in other tumor datasets. In TCGA tumors examined by RNA-Seq (Additional File 1 Figure S6b) as in other tumors examined by microarrays (Additional File 1, Figure S6c), high expression of the 8 DEGs reduced OS specifically in TNBCs. Altogether, this suggests that the high ERRα-ZEB1 activity (proxied by expression of the 8 DEGs) constitutes a factor of poor prognosis specifically in TNBC.

Transcription factors act in cooperation with co-regulators (co-activators or co-repressors) that can impact on their target repertoire and thereby influence their eventual physiological activities. For instance, the ERRα nuclear receptor is highly involved in the regulation of cellular metabolism in tissues with high energy demands such as liver or muscle [29]. This is accomplished via binding of the receptor to genomic regions of metabolic genes and recruitment of members of the PGC-1 (α and/or β) family of co-activators [5, 24]. However, ERRα does not regulate these metabolic genes in breast cancer cell lines although ChIP-Seq experiments demonstrate its presence on their promoters [6, 7]. In contrast, ERRα directly activates the expression of genes involved in parameters of cancer progression in a PGC-1-independent manner, thus questioning the identity of the co-regulators involved in these regulations. Various methods have been used to identify such ERRα co-factors. Some, such as LSD1, have been identified through a candidate gene approach [12], others, such as RIP140 based on their physical interaction with ERRα [17]. Although a recent publication documents the identification of the ERRα interactome in mouse liver cells [30], there is a general lack of unbiased methods to identify ERRα co-regulators in non-metabolic tissues.

We here applied a bio-informatics approach, sPLS, to the analysis of publically available RNA-Seq data from breast cancers, aiming to identify factors that could contribute to the expression of ERRα direct target genes. These were first identified in a breast cancer cell line by an RNA-Seq approach (this report) that we compared to our published ERRα ChIP-Seq [7]. We further focused on the genes appearing in the principal GO term, cell migration, a parameter that is regulated by ERRα [6, 10, 12, 18, 31]. The sPLS method has been already validated on the study of RNA-Seq from breast cancer cell lines [7]. Its reliability was here insured by performing the analysis three times each with ten replicates, producing similar results. Our analysis suggested ZEB1 as a main potential ERRα transcriptional co-regulator in breast cancers, including in TNBC. Interestingly, further analysis showed that ZEB1 appeared in first position in all TNBC subtypes defined by Jiang et al. [19]. Together with ERRα, this factor is suggested to regulate the expression of 8 out of the 45 DEGs that were used for the sPLS analysis, in nearly all tumor types. This hypothesis was experimentally confirmed in various TNBC cell lines. Interestingly, this effect is ZEB1-specific, since depletion of the ERRα co-activator PGC-1α did not influence the expression of the 8 DEGs. In addition, further experiments indicated that some, but not all, other ZEB1 targets could also be regulated by ERRα in various TNBC cell lines. Altogether, this suggests that, in TNBC, ERRα and ZEB1 co-regulate a limited number of targets, which may however extend above 8.

ZEB1 displays different roles in cancer [27]. One of its main reported functions is to induce EMT, which is accomplished via various molecular activities including repression of its target genes. This is performed in a direct manner via recruitment on ZEB1-response element (E-boxes) or an indirect manner via the reciprocal regulation of the expression of miR-200 family members, which in turn impact key actors of EMT [25, 26, 32]. Strikingly, ZEB1 can also activate the expression of its target genes. To this end, ZEB1 physically interacts with various transcription factors, indirectly binds to their respective response elements that are distinct from E-boxes and promotes transcriptional activation of a limited subset of these target genes. ZEB1 co-activating activities have been documented to occur with YAP1, LEF1 and JUN-FOSL1, depending on the cell type considered [23, 33–34]. We here demonstrate a similar situation with ZEB1 interacting with ERRα, indirectly binding to ERREs and stimulating the expression of a subset of ERRα target genes. Altogether, this shows that, in addition to its autonomous activities, ZEB1 acts as a co-regulator that potentiates the effects of various transcription factors on a limited amount of targets. Interestingly, all these heteronomous activities converge to promote cancer progression. In the present case, this is reflected by the reported cancer functions of the ZEB1-ERRα targets documented here. Indeed, ADAMTS12 is a metalloprotease that favors cell migration in breast cancer cells [35]; CDH5, a transmembrane protein, is a marker for metastatic breast cancers [36]; DDR2 is a receptor tyrosine kinase that stabilizes EMT [37]; ETS1 is a transcription factor that acts in TNBC progression [38]; MCU, a calcium transporter, increases angiogenesis in metastatic niches [39]; PECAM1 is a cell adhesion molecule that induces EMT [40]; PLXNC1, a receptor for semaphorins is involved in breast cancer progression [41]; SPARC is an ECM protein that promotes lung metastasis of breast cancers [42]. These activities are also consistent with the various cancer promoting effects of ERRα [18, 43]. Interestingly, our data show that expression of all DEGs (except MCU) is linked to the EMT score of TCGA breast cancers.

Identify markers that would predict the survival of breast cancer patients, particularly in the case of TNBC, the most aggressive breast cancer type, is a key issue. Indeed this could influence the choice of a therapy type used and thereby possibly minimize the unwanted negative effects of these treatments. On another hand, if the pro-cancer activities of these markers are reported, they could be used as targets, provided that they are expressed in the considered cases, reinforcing the need to study their expression.

High ZEB1 expression is associated with aggressiveness, metastasis and resistance to therapy [27, 44]. As such, it has been associated to a poor prognosis in pancreatic, ovarian and breast cancer [23, 45–46]. In contrast, our analysis indicates that high ZEB1 expression is linked to overall survival through a non-significant trend and selectively in TNBC. The reason for this discrepancy is unknown. However, it is possible that the tumors used as a whole for RNA extraction and studied for gene expression have a different content in non-tumor cells (e.g. immune or stromal cells). In turn, these different distributions may influence the results of the transcriptome analysis and thus bias the correlation to survival. In this line, it should be noted that ZEB1 expression is restricted to a subpopulation of stem-like malignant cells in pancreatic cancers [45]. On another hand, high ZEB1 expression in tumor-associated macrophages of ovarian carcinomas constitutes a factor of poor prognosis [46]. However, it should be noted that the published analysis performed in breast cancer considered distant metastasis free survival performed on microarray-based transcriptome [23]. In contrast, we here analyze overall survival examined on RNA-Seq-based expression data. Overall, this also suggests that the activity of ZEB1 (and not only the expression of its corresponding RNAs) should be studied and tentatively correlated with survival. This reasoning can also be applied to ERRα. Indeed, the expression of its mRNA is not different between tumors and normal breast tissue [47]. Furthermore, a recent publication shows that examining the ERRα transcriptional target repertoire in thyroid tumors allows to discriminate between two tumor subtypes [48]. Study of the ERRα-ZEB1 DEGs shows this hypothesis to be conclusive.

Mortality in breast cancer patients of African American patients is higher than that observed in White women [49]. This phenomenon can be due to socioeconomic- as well as biological factors [50]. Moreover, TNBC are overrepresented in African American women as compared to White patients. This is also the case for pre-menopausal vs post-menopausal patients [51], with pre-menopausal cancers being more aggressive than post-menopausal ones. There is thus a high interest in identifying marker genes of sensitivity/survival in African Americans or pre-menopausal women [52–53]. Here we find that high expression of several ERRα-ZEB1 DEGs can be explained by ZEB1 expression in pre-menopausal patients and, in a more moderated manner, in African American women. Strikingly, this applies to TNBC patients, suggesting that expression of the 8 DEGs could be used for survival prediction in these groups although the low number of TNBC sorted according to these parameters is too low to allow a survival analysis.

In summary, we here show that the ERRα-ZEB1 complex is a major actor in breast cancer progression in particular for TNBC. Furthermore, high expression of the 8 ERRα-ZEB1 DEGs, when considered all together, constitutes a prognosis factor (in terms of overall survival) for these highly aggressive tumors.

Transcription Factor

ERRα

Estrogen-Related Receptorα

TNBC

Triple-negative Breast Cancer

Glucocorticoid Receptor

Transcriptional Regulator

sPLS

sparse Partial Least Square

DEG

Differentially Expressed Gene

EMT

Epithelium-to-Mesenchyme Transition

Gene Ontology

ERRE

ERR-Response Element

Overall Survival

Ethics approval and consent to participate.

Not applicable

Consent for publication.

Not applicable

Data and material availability.

The datasets generated during and/or analyzed during the current study as well as all materials are available from the corresponding author on reasonable request.

Competing Interests

OT received honoraria for lectures, presentations, and advisory board from Roche, Pfizer, Novartis-Sandoz, Lilly, MSD, Astra-Zeneca, Pierre Fabre, Seagen, Daiichi-Sankyo, Gilead, Eisai, Menarini-Stemline, Veracyte and Exact Sciences. Funders had no role in the design of the study, collection and analysis of data and decision to publish.

Funding.

This work was funded by Ligue contre le Cancer (comité Ardèche and Drôme) and JoRiss/ENS research programs. JRS is funded by the Chinese Scholarship Council (CSC).

Author Contribution

JRS and JMV designed the experiments, JRS, CP and CC performed experiments, CC, OT, MLR, TLS and JMV analyzed data, JRS and JMV wrote the paper.

Acknowledgement

We are indebted to the Plateforme de Séquençage IGFL (PSI), in particular to Sandrine Hughes and Benjamin Gillet, for the sequencing of the RNA-Seq.

Data Availability

Results of the RNA-Seq are available at GEO website under accession number GSE198923

Dey A, Varelas X, Guan KL. Targeting the Hippo pathway in cancer, fibrosis, wound healing and regenerative medicine. Nat Rev Drug Discov. 2020;19(7):480–94.
Koontz LM, Liu-Chittenden Y, Yin F, Zheng YG, Yu JZ, Huang B, et al. The Hippo effector Yorkie controls normal tissue growth by antagonizing scalloped-mediated default repression. Dev Cell. 2013;25(4):388–401.
Chinenov Y, Gupte R, Rogatsky I. Nuclear receptors in inflammation control: repression by GR and beyond. Mol Cell Endocrinol. 2013;380(1–2):55–64.
Giguère V. Transcriptional control of energy homeostasis by the estrogen-related receptors. Endocr Rev. 2008;29(6):677–96.
Villena JA, Kralli A. ERRalpha: a metabolic function for the oldest orphan. Trends Endocrinol Metab. 2008;19(8):269–76.
Sailland J, Tribollet V, Forcet C, Billon C, Barenton B, Carnesecchi J, et al. Estrogen-related receptor α decreases RHOA stability to induce orientated cell migration. Proc Natl Acad Sci U S A. 2014;111(42):15108–13.
Cerutti C, Zhang L, Tribollet V, Shi JR, Brillet R, Gillet B, et al. Computational identification of new potential transcriptional partners of ERRα in breast cancer cells: specific partners for specific targets. Sci Rep. 2022;12(1):3826.
Cerutti C, Shi JR, Vanacker JM. Multifaceted Transcriptional Network of Estrogen-Related Receptor Alpha in Health and Disease. Int J Mol Sci. 2023;24(5):4265.
Bianco S, Lanvin O, Tribollet V, Macari C, North S, Vanacker JM. Modulating estrogen receptor-related receptor-alpha activity inhibits cell proliferation. J Biol Chem. 2009;284(35):23286–92.
Dwyer MA, Joseph JD, Wade HE, Eaton ML, Kunder RS, Kazmin D, et al. WNT11 expression is induced by estrogen-related receptor alpha and beta-catenin and acts in an autocrine manner to increase cancer cell migration. Cancer Res. 2010;70(22):9298–308.
Zou C, Yu S, Xu Z, Wu D, Ng CF, Yao X, et al. ERRα augments HIF-1 signalling by directly interacting with HIF-1α in normoxic and hypoxic prostate cancer cells. J Pathol. 2014;233(1):61–73.
Carnesecchi J, Forcet C, Zhang L, Tribollet V, Barenton B, Boudra R, et al. ERRα induces H3K9 demethylation by LSD1 to promote cell invasion. Proc Natl Acad Sci U S A. 2017;114(15):3909–14.
Tribollet V, Cerutti C, Géloën A, Berger E, De Mets R, Balland M, et al. ERRα coordinates actin and focal adhesion dynamics. Cancer Gene Ther. 2022;29(10):1429–38.
Villena JA, Hock MB, Chang WY, Barcas JE, Giguère V, Kralli A. Orphan nuclear receptor estrogen-related receptor alpha is essential for adaptive thermogenesis. Proc Natl Acad Sci U S A. 2007;104(4):1418–23.
Heck S, Rom J, Thewes V, Becker N, Blume B, Sinn HP, et al. Estrogen-related receptor alpha expression and function is associated with the transcriptional coregulator AIB1 in breast carcinoma. Cancer Res. 2009;69(12):5186–93.
Liu D, Benlhabib H, Mendelson CR. cAMP enhances estrogen-related receptor alpha (ERRalpha) transcriptional activity at the SP-A promoter by increasing its interaction with protein kinase A and steroid receptor coactivator 2 (SRC-2). Mol Endocrinol. 2009;23(6):772–83.
Castet A, Herledan A, Bonnet S, Jalaguier S, Vanacker JM, Cavaillès V. Receptor-interacting protein 140 differentially regulates estrogen receptor-related receptor transactivation depending on target genes. Mol Endocrinol. 2006;20(5):1035–47.
Vanacker JM, Forcet C. ERRα: unraveling its role as a key player in cell migration. Oncogene. 2024;43(6):379–87.
Jiang YZ, Ma D, Suo C, Shi J, Xue M, Hu X, et al. Genomic and Transcriptomic Landscape of Triple-Negative Breast Cancers: Subtypes and Treatment Strategies. Cancer Cell. 2019;35(3):428–e4405.
Gambardella G, Viscido G, Tumaini B, Isacchi A, Bosotti R, di Bernardo D. A single-cell analysis of breast cancer cell lines to study tumour heterogeneity and drug response. Nat Commun. 2022;13(1):1714.
Parker JS, Mullins M, Cheang MC, Leung S, Voduc D, Vickery T, et al. Supervised risk predictor of breast cancer based on intrinsic subtypes. J Clin Oncol. 2009;27(8):1160–7.
Prat A, Parker JS, Karginova O, Fan C, Livasy C, Herschkowitz JI, et al. Phenotypic and molecular characterization of the claudin-low intrinsic subtype of breast cancer. Breast Cancer Res. 2010;12(5):R68.
Feldker N, Ferrazzi F, Schuhwerk H, Widholz SA, Guenther K, Frisch I, et al. Genome-wide cooperation of EMT transcription factor ZEB1 with YAP and AP-1 in breast cancer. EMBO J. 2020;39(17):e103209.
Deblois G, St-Pierre J, Giguère V. The PGC-1/ERR signaling axis in cancer. Oncogene. 2013;32(30):3483–90.
Burk U, Schubert J, Wellner U, Schmalhofer O, Vincan E, Spaderna S, et al. A reciprocal repression between ZEB1 and members of the miR-200 family promotes EMT and invasion in cancer cells. EMBO Rep. 2008;9(6):582–9.
Gregory PA, Bert AG, Paterson EL, Barry SC, Tsykin A, Farshid G, et al. The miR-200 family and miR-205 regulate epithelial to mesenchymal transition by targeting ZEB1 and SIP1. Nat Cell Biol. 2008;10(5):593–601.
Caramel J, Ligier M, Puisieux A. Pleiotropic Roles for ZEB1 in Cancer. Cancer Res. 2018;78(1):30–5.
Polyak K, Weinberg RA. Transitions between epithelial and mesenchymal states: acquisition of malignant and stem cell traits. Nat Rev Cancer. 2009;9(4):265–73.
Xia H, Dufour CR, Giguère V. ERRalpha as a Bridge Between Transcription and Function: Role in Liver Metabolism and Disease. Front Endocrinol (Lausanne). 2019;10:206.
Scholtes C, Dufour CR, Pleynet E, Kamyabiazar S, Hutton P, Baby R, et al. Identification of a chromatin-bound ERRα interactome network in mouse liver. Mol Metab. 2024;83:101925.
Zhang L, Carnesecchi J, Cerutti C, Tribollet V, Périan S, Forcet C, et al. LSD1-ERRα complex requires NRF1 to positively regulate transcription and cell invasion. Sci Rep. 2018;8(1):10041.
Brabletz S, Bajdak K, Meidhof S, Burk U, Niedermann G, Firat E, et al. The ZEB1/miR-200 feedback loop controls Notch signalling in cancer cells. EMBO J. 2011;30(4):770–82.
Lehmann W, Mossmann D, Kleemann J, Mock K, Meisinger C, Brummer T, et al. ZEB1 turns into a transcriptional activator by interacting with YAP1 in aggressive cancer types. Nat Commun. 2016;7:10498.
Rosmaninho P, Mükusch S, Piscopo V, Teixeira V, Raposo AA, Warta R, et al. Zeb1 potentiates genome-wide gene transcription with Lef1 to promote glioblastoma cell invasion. EMBO J. 2018;37(15):e97115.
Mohamedi Y, Fontanil T, Cal S, Cobo T, Obaya ÁJ. ADAMTS-12: Functions and Challenges for a Complex Metalloprotease. Front Mol Biosci. 2021;8:686763.
Fry SA, Robertson CE, Swann R, Dwek MV. Cadherin-5: a biomarker for metastatic breast cancer with optimum efficacy in oestrogen receptor-positive breast cancers with vascular invasion. Br J Cancer. 2016;114(9):1019–26.
Trono P, Ottavi F. Rosano' L. Novel insights into the role of Discoidin domain receptor 2 (DDR2) in cancer progression: a new avenue of therapeutic intervention. Matrix Biol. 2024;12:31–9.
Dittmer J. The role of the transcription factor Ets1 in carcinoma. Semin Cancer Biol. 2015;35:20–38.
Zhang YY, Kong LQ, Zhu XD, Cai H, Wang CH, Shi WK, et al. CD31 regulates metastasis by inducing epithelial-mesenchymal transition in hepatocellular carcinoma via the ITGB1-FAK-Akt signaling pathway. Cancer Lett. 2018;429:29–40.
Zheng X, Lu S, He Z, Huang H, Yao Z, Miao Y, et al. MCU-dependent negative sorting of miR-4488 to extracellular vesicles enhances angiogenesis and promotes breast cancer metastatic colonization. Oncogene. 2020;39(46):6975–89.
Aiyappa-Maudsley R, McLoughlin LFV, Hughes TA. Semaphorins and Their Roles in Breast Cancer: Implications for Therapy Resistance. Int J Mol Sci. 2023;24(17):13093.
Lepucki A, Orlinska K, Mielczarek-Palacz A, Kabut J, Olczyk P, Komosinska-Vassev K. The Role of Extracellular Matrix Proteins in Breast Cancer. J Clin Med. 2022;11(5):1250.
Bianco S, Sailland J, Vanacker JM. ERRs and cancers: effects on metabolism and on proliferation and migration capacities. J Steroid Biochem Mol Biol. 2012;130(3–5):180–5.
Stemmler MP, Eccles RL, Brabletz S, Brabletz T. Non-redundant functions of EMT transcription factors. Nat Cell Biol. 2019;21(1):102–12.
Wellner U, Schubert J, Burk UC, Schmalhofer O, Zhu F, Sonntag A, et al. The EMT-activator ZEB1 promotes tumorigenicity by repressing stemness-inhibiting microRNAs. Nat Cell Biol. 2009;11(12):1487–95.
Cortés M, Sanchez-Moral L, de Barrios O, Fernandez-Acenero MJ, Martinez-Campanario MC, Esteve-Codina A, et al. Tumor-associated macrophages (TAMs) depend on ZEB1 for their cancer-promoting roles. EMBO J. 2017;36(22):3336–55.
Carnesecchi J, Cerutti C, Vanacker JM, Forcet C. ERRα protein is stabilized by LSD1 in a demethylation-independent manner. PLoS ONE. 2017;12(11):e0188871.
Chen W, Song YS, Lee HS, Lin CW, Lee J, Kang YE, et al. Estrogen-related receptor alpha promotes thyroid tumor cell survival via a tumor subtype-specific regulation of target gene networks. Oncogene. 2024;43(31):2431–46.
Newman LA, Griffith KA, Jatoi I, Simon MS, Crowe JP, Colditz GA. Meta-analysis of survival in African American and white American patients with breast cancer: ethnicity compared with socioeconomic status. J Clin Oncol. 2006;24(9):1342–9.
Dietze EC, Sistrunk C, Miranda-Carboni G, O'Regan R, Seewaldt VL. Triple-negative breast cancer in African-American women: disparities versus biology. Nat Rev Cancer. 2015;15(4):248–54.
Zagami P, Carey LA. Triple negative breast cancer: Pitfalls and progress. NPJ Breast Cancer. 2022;8(1):95.
Martini R, Delpe P, Chu TR, Arora K, Lord B, Verma A, et al. African Ancestry-Associated Gene Expression Profiles in Triple-Negative Breast Cancer Underlie Altered Tumor Biology and Clinical Outcome in Women of African Descent. Cancer Discov. 2022;12(11):2530–51.
Kumar B, Khatpe AS, Guanglong J, Batic K, Bhat-Nakshatri P, Granadir MM, et al. Stromal heterogeneity may explain increased incidence of metaplastic breast cancer in women of African descent. Nat Commun. 2023;14(1):5683.
Poulard C, Jacquemetton J, Pham TH, Le Romancer M. Using proximity ligation assay to detect protein arginine methylation. Methods. 2020;175:66–71.
Davis SR, Baber RJ. Treating menopause - MHT and beyond. Nat Rev Endocrinol. 2022;18(8):490–502.
Chakraborty P, George JT, Tripathi S, Levine H, Jolly MK. Comparative Study of Transcriptomics-Based Scoring Metrics for the Epithelial-Hybrid-Mesenchymal Spectrum. Front Bioeng Biotechnol. 2020;8:220.
Győrffy B. Integrated analysis of public datasets for the discovery and validation of survival-associated genes in solid tumors. Innov (Camb). 2024;5(3):100625.
Lehmann BD, Bauer JA, Chen X, Sanders ME, Chakravarthy AB, Shyr Y, et al. Identification of human triple-negative breast cancer subtypes and preclinical models for selection of targeted therapies. J Clin Invest. 2011;121(7):2750–67.

Competing interest reported. OT received honoraria for lectures, presentations, and advisory board from Roche, Pfizer, Novartis-Sandoz, Lilly, MSD, Astra-Zeneca, Pierre Fabre, Seagen, Daiichi-Sankyo, Gilead, Eisai, Menarini-Stemline, Veracyte and Exact Sciences. Funders had no role in the design of the study, collection and analysis of data and decision to publish.

Download PDF

Version 1

posted

You are reading this latest preprint version

An ERRα-ZEB1 transcriptional signature predicts survival in triple-negative breast cancers

Status:

Version 1

Abstract

Background.

Methods.

Results.

Conclusions.

Figures

Introduction

Materials and Methods

Results

Discussion

Abbreviations

Declarations

Data and material availability.

Competing Interests

Funding.

Author Contribution

Acknowledgement

Data Availability

References

Additional Declarations

Supplementary Files

Status:

Version 1