Adverse Crosstalk Between Extracellular Matrix Remodeling and Ferroptosis in Basal Breast Cancer

1. Introduction

The breast cancer is the most common cancer newly diagnosed for women in United States in 2020 [1]. Breast cancer is a heterogeneous disease with different molecular subtypes defined by distinct molecular classes associated to prognosis: claudin-low, normal-like, luminal A, luminal B, HER2 and basal [2] and confirmed by gene quantification in pam50 classification [3]. Breast cancer is the second leading cause of mortality in woman due to rapid metastasis and disease recurrence [4]. The breast tissue is in a unique microenvironment with plentiful adipocytes infiltrating. Previous studies have shown that adipocytes can regulate fatty acid metabolism, enhance the invasion and metastasis of breast cancer [5]. Ferroptosis is an iron-dependent regulated form of cell death caused by the accumulation of lipid-based reactive oxygen species (ROS) [6]. Prerequisites for ferroptosis include iron metabolism, mitochondrial metabolism, synthesis of polyunsaturated fatty acid phospholipid (PUFA-PL) and lipid peroxidation [7]. Targeting ferroptosis have been proposed to treat breast cancer. There is increasing evidence that ferroptosis can inhibit breast cancer cell growth, improve the sensitivity of chemotherapy and radiotherapy and inhibit distant metastases [5]. From 2001 to 2003, the Stockwell Lab performed a screen to identify compounds that kill cells engineered to be tumourigenic (harbouring the RAS mutant), without killing their isogenic parental precursors. One of the most efficient compounds was identified and named “erastin” after its ability to “Eradicate RAS-and Small T transformed cells” [8]. Subsequently, they identified RSL3, which was also named after its “oncogenic-RAS-selective lethal” property in 2008 [9]. Small molecule–induced ferroptosis has been shown to have a strong inhibitory effect on tumor growth in a drug-resistant environment, which may increase the sensitivity of the tumor to chemotherapeutic treatment [10]. Ferroptosis is also considered an important cell death mechanism caused by a number of therapies, including chemotherapy, radiotherapy (RT), targeted therapy and immunotherapy [11], but in opposite way, the tumor cells with ferroptosis could diminish anti-tumor immune response by inhibiting the antigen presenting cells [12]. The tumor microenvironment (TME) plays a notable role in cancer progression. It includes pH and oxygen levels, the extracellular matrix (ECM), connective tissue, infiltrating immune cells, and the vasculature of the tumor. Interaction between the ECM and the tumor cells activates key signaling pathways that promote tumor proliferation, invasion, and metastasis. This notably influences many tumors as the ECM can comprise up to 60% of the tumor mass [13].

In the present work, by a text-mining approach integrated in transcriptome experiments, a link between Ferroptosis and ECM remodeling was done through gene related regulation in adverse prognosis of breast cancer but also in TNBC cellular model stimulated by ferroptosis activators.

2. Materials and Methods

• Determination of ferroptosis in breast related genes

Using keyword “Ferroptosis in breast” keyword, a co-occurrence of citations with coding gene identifiers was searched in the article abstracts of the PUBMED database was searched with the "Génie" algorithm [14]. Bioinformatics analyses were realized in R software environment version 4.2.1. False Discovery Rate (FDR) of the significant associated genes were analyzed by qqplot visualization with qqman R package version 0.1.8 [15]. Further investigations of text-mining associations with gene identifiers have been confirmed with the “GeneValorization” application [16] on the National Center for Biotechnolgy Information (NCBI) database [17]. The results of this text-mining were drawn as circleplot of gene-keywords co-occurrence with circlize R-package version 0.4.15 and as alluvial plot.

• Transcriptome cohort of breast cancer for patients treated with anthracyclin and taxanes

Transcriptome normalized matrix of dataset GSE25066 [18] was downloaded at the following address: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=gse25066 (accessed on 2023, may 18^th) and annotated with the corresponding technology platform GPL96 [HG-U133A] Affymetrix Human Genome U133A Array available at the following address: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL96 (accessed on 2023, may 18^th).

• Transcriptome dataset testing effect of ferroptosis activators on MDA-MB-231 triple negative breast cancer cell line

Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced (FPKM) Transcript quantification performed by original pipeline of dataset GSE173905 [19] were downloaded at the following address:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE173905 (accessed on 2023, may 18^th). After sequencing on Illumina NovaSeq 6000 technology, the original pipeline aligned reads on human genome with the reference genome was built using Hisat2 v2.0.5 [20] and paired-end clean reads were aligned to the reference genome using Hisat2 v2.0.5 software and FPKM, expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced were computed on counts obtained with Feature Counts v1.5.0-p3 software [21]. On selected genes for Ferroptosis/Extracellular matrix remodeling signature, unsupervised principal component analysis was performed with FactoMiner R package version 2.8 [22].

• Gene expression analyses and association to the breast cancer prognosis

Distant relapse free survival (DRFS) from dataset GSE25066 [18] was used as outcome to performed iterative univariate Cox analysis against expression of the genes identified as being significantly associated with the keyword "ferroptosis in breast". During trial follow-up of adjuvant therapy in breast cancer, distant relapse–free survival (DRFS) could be used as endpoint [18]. DRFS was defined as the interval from initial diagnostic biopsy until diagnosis of distant metastasis or death from breast cancer, non-breast cancer, or unknown causes [23]. This iteration of univariate Cox analysis was automatized with loopcolcox R-package version 1.0.0 available at the following address: https://github.com/cdesterke/loopcolcox (accessed on 2023, may 18^th). Univariate Kaplan Meier and survival optimal threshold on variables were performed with survminer R-package version 0.4.9 and survival R-package version 3.3.1. On genes with adverse prognosis association a functional enrichment was performed with CPG signature from MsigDb database [24] through Toppgene online application [25]. A breast cancer related signature network was drawn with Cytoscape standalone software version 3.9.1 [26]. An expression molecular score related to Ferroptosis/Extracellular matrix remodeling functionalities was done by computing the sum of the product between Cox beta coefficients and expression of the eleven selected genes. For the eleven genes belonging to the Ferroptosis/Extracellular matrix remodeling signature a multi-ROC analysis was done against estrogen/progesterone receptors status detected in immunohistochemistry with the R-package multirocauc version 1.0.0 available at the address: https://github.com/cdesterke/multirocauc (accessed on 2023, may 18^th) (Supplemental Figure 1A). A multivariate Cox model was built with DRFS as outcome and with incorporation of molecular score and relevant clinical parameters. This DRFS multivariate model was assessed by testing linearity of residuals at global level and for each individual included parameters with Schoenfeld tests. Calibration of the DRFS multivariate model at 10 months of follow-up was done by 500 iterations of bootstrap with rms R-package version 6.7.0. Nomogram validated at 10 months of follow-up was drawn for DRFS multivariate model with regplot R-package version 1.1.

3. Results

3.1. Ferroptosis gene expression associated to the prognosis of patients with breast cancer

A text-mining approach was employed to found genes related to the ferroptosis in breast literature. The text-mining algorithm “Génie” was employed by querying Pubmed with “ferroptosis in breast” keywords. This query returns a list of 910 individual genes with significant False Discovery Rate (FDR) and positive in at least 10 distinct articles (Figure 1A). Iterations of univariate Cox analyses against distant relapse free survival (DRFS) outcome of patients from transcriptome dataset GSE25066 [18] was done for each of the 910 genes selected by text-mining. For best fifty ranked genes associated in their expression to the prognosis of the patients, some of them were found favorable and other adverse according the sense of their beta coefficients or hazard ratios (Figure 1B). For example, PTTG1 (PTTG1 regulator of sister chromatid separation, securin) is the ferroptosis in breast related gene found with the most adverse expression in breast cancer samples (hazard ratio=2.21) followed by ENO1 (enolase 1) with hazard ratio 2.23, SLC7A5 (solute carrier family 7 member 5) with hazard ratio 1.51, ITCH (itchy E3 ubiquitin protein ligase) with hazard ratio 3.86, EGLN1 (egl-9 family hypoxia inducible factor 1) with hazard ratio 2.26 (Figure 1B and Supplemental Table 1).

Among the 252 significant genes associated to the prognosis, a filtration on positivity of Cox beta-coefficients was done to retain 171 genes associated to adverse prognosis. Functional enrichment performed with these 171 adverse genes on MsigDb CPG signature database highlighted main enrichment with of these genes in breast cancer published signatures such as: SMID-breast cancer basal up [27] and SOTIRIOU-breast cancer grade 1 vs 3 up [28] (Figure 2A). A breast cancer network was drawn with the ferroptosis. Some ferroptosis-related genes were found shared between the two independent breast cancer transcriptome cohort signatures (Figure 2B), such as: FOXM1 (forkhead box M1), SLC7A5 (solute carrier family 7 member 5), CHEK1 (checkpoint kinase 1), MKI67 (marker of proliferation Ki-67), TTK (TTK protein kinase), PTTG1 (PTTG1 regulator of sister chromatid separation), BIRC5 (baculoviral IAP repeat containing 5), MELK (maternal embryonic leucine zipper kinase), CCNA2 (cyclin A2), NEK2 (NIMA related kinase 2), TPX2 (TPX2 microtubule nucleation factor), CDK1 (cyclin dependent kinase 1) and CDC20 (cell division cycle 20). These results suggest that the ferroptosis in breast related text-mining approach allowed to predict breast cancer transcriptome experiments already published and among these genes some are associated to adverse prognosis of the patients. To validate Génie text-mining approach an independent text-mining algorithm “GeneValorization” was employed still to query PUBMED database with “ferroptosis” keyword but also some others with relevance in the context of the study such as: cancer stem cell (CSC), Extracellular matrix (ECM) remodeling, breast cancer, lipid peroxidation and regulated cell death. This validation allowed to highlight best 15 top ranked genes in text-mining sharing these keywords associations (Figure 2C): SET (S (ET nuclear proto-oncogene), TNF (tumor necrosis factor), HMOX1 (heme oxygenase 1), IL6 (interleukin 6), TRFC (transferrin receptor), ATF4 (activating transcription factor 4), HMGB1 (high mobility group box 1), KRAS (KRAS proto-oncogene, GTPase), EGFR (epidermal growth factor receptor), TLR4 (toll like receptor 4), HIF1A (hypoxia inducible factor 1 subunit alpha), ATG5 (autophagy-related 5), LCN2 (lipocalin 2), CDKN2A (cyclin dependent kinase inhibitor 2A), MET (MET proto-oncogene, receptor tyrosine kinase).

3.2. Ferroptosis/extracellular matrix remodeling expression score and prognosis in breast cancer

Among ferroptosis validated genes in breast cancer (Figure 2C), eleven of them harbored some association with extracellular matrix remodeling literature in PUBMED database, especially: TNF, SET and IL6 were found highly cited in this specific literature (Figure 3A). Exploring individual expression of these eleven genes against immunohistochemistry status of estrogen receptor highlighted great importance of expression for CDKN2A (aera under curve (AUC: 0.76) and EGFR (AUC: 0.79) (Figure 3B) for this hormonal prediction. Unsupervised principal component analysis based on expression of the eleven ferroptosis/extracellular matrix remodeling-related genes well stratified patient groups according their estrogen receptor status on first principal axis (Supplemental Figure 2A). Concerning immunohistochemistry status of progesterone receptor, expression of CDKN2A appeared decisive (AUC: 0.72, Figure 3C).

Unsupervised principal component analysis based on expression of the eleven ferroptosis/extracellular matrix remodeling-related genes well stratified groups of patients according their progesterone receptor status on first principal axis (Supplemental Figure 2B). Similarly, principal component analysis gradually stratified tumor grades on the first principal axis (Figure 4A), and pam50 molecular groups (Figure 4B) suggesting that the eleven gene signature related to ferroptosis/ECM remodeling could be interesting for prognostication of the breast cancer patients. A molecular score was computed on the expression of the eleven ferroptosis/extracellular matrix remodeling related genes. Optimal cutpoint was determined on DRFS residuals of the molecular score (Figure 4C). This threshold cutpoint was found at 48.74 to stratify the breast cancer cohort in two groups. No significant difference was found between the two groups of patients harboring low and high values of ferroptosis/ECM remodeling molecular score (p-value=0.52, Table 1). Concerning immunohistochemistry status of estrogen receptor, a significant higher proportion of negative patients was found in group of patients with high values of molecular score (p-value<1E-4, Table 1), and it was the same observation for progesterone receptor status (p-value<1E-4, Table 1). Concerning pam50 molecular classification, groups of patients with high values of ferroptosis/ECM remodeling score present higher proportion of basal type samples (p-value<1E-4, Table 1). No significant difference was observed on tumor stages between the two groups of patients (p-value=0.38, Table 1) but a higher proportion of patients N3 for nodal status was observed in group of patients with high value of molecular score (p-value=0.013, Table 1). Concerning clinical AJCC staging a significant difference was observed with increase proportion of stages IIIA, IIIB, IIIC and inflammatory in group of patients with high value of molecular score (p-value=0.03, Table 1). As observed by PCA on the expression of the eleven genes which composed the molecular score (Figure 4A), the grade variable was found significant between the two groups of patients (p-value<1E-4, Table 1). DRFS status and time parameters were also confirmed as significant between the two groups of patients (p-value<1E-4, Table 1). Effectively, Kaplan Meier with DRFS censor stratified on ferroptosis/ECM remodeling is very significant (Figure 4D), with a major worst prognosis for patients which harbored a molecular score over 48.74 as threshold. Type of taxanes (Taxol, Taxotere) administrated during follow-up of the patients present no associations with the Ferroptosis/ECM remodeling groups of patients (p-value=0.69, Table 1).

Relevant clinical parameters were integrated in a multivariable Cox model censored on DRFS with group stratification based on ferroptosis/ECM remodeling molecular score. This multivariable model which harbored a concordant index of 0.77 was found highly significant by likelihood ratio test (p-value=2E-14) (Figure 5A). Global and individual Schoenfeld test attested linear distribution residuals from included parameters: age of patients, nodal status, pam50 molecular classification, grading, and ferroptosis/ECM remodeling molecular score (Supplemental Figure 3). In this multivariable model high values of nodal status (N1 and N2,3) were found as independent adverse parameter of DRFS (N1 versus N0 hazard ratio: 2.20, p-value=5.91E-3, N23 versus N0 hazard ratio: 3.39, p-value=7.15E-5, Table 2 and Figure 5A). Among molecular classification, basal subtype appeared as adverse group of breast cancer with normal like as reference (Hazard ratio:2.92, p-value 4.45E-2, Table 2 and Figure 5A). In the DRFS multivariate model, ferroptosis/ECM remodeling molecular score appeared as an adverse independent parameter in the prognosis of breast cancer patients (high score vs low score hazard ratio: 2.69, p-value=1.17E-5). The multivariate model could be calibrated at ten months of follow-up with five hundred iterations by Kaplan Meier method (Figure 5B): this calibration showed that the multivariable model is stable at 10 months of follow-up because near of ideal model in gray. The corresponding nomogram of the model was drawn for a prediction at 10 months of follow up (Figure 5C). This representation confirms the important part of the molecular score in the multivariate model. Indeed, at 10 months of follow-up, the molecular score appears dispersed between the range of point values of the model (10-70) for a DRFS probability at 10 months comprised between (0.006 and 0.1) (Figure 5C).

3.3. Ferroptosis/ECM remodeling signature is regulated by ferroptosis activators in triple negative breast cancer cells

To verify ferroptosis functional of gene members belonging to the molecular score, transcriptome of MDA-MB-231 cells (cellular model of Triple Negative Breast Cancer (TNBC)) stimulated by two distinct ferroptosis activators: erastin and RSL3 was analyzed (GSE173905) [19]. Based on the expression of the eleven genes belonging to ferroptosis/ECM molecular score, an unsupervised principal component analysis was performed with samples of GSE173905. This multivariate analysis confirmed that gene members composing the molecular score are regulated by the two distinct ferroptosis activators: erastin and RSL3 (p-value=1.08E-5, Figure 6A). During this ferroptosis functional activations, gene members are regulated by distinct patterns (Figure 6B). Some genes like HIF1A (Figure 6C) and LCN2, EGFR, IL6 (Figure 6D) were found up regulated by both ferroptosis activators. Some genes like HMGB1, KRAS, SET (Figure 6E) were found down regulated by both ferroptosis activators. MET and TLR4 were found regulated in opposite way by the two ferroptosis activators (Figure 6F). Finally, TNF and CDKN2A were not found detected by sequencing in this MDA-MB-231 TNBC cellular model.

4. Discussion

During this study, by text mining approach, a list of genes was defined as related to ferroptosis cellular functionality known as important way cellular death implicated in tumor response to therapies [10]. Surprisingly, in transcriptome of breast tumors under therapies (GSE173905) [18], the majority of ferroptosis-related genes were found having expression associated with adverse distant relapse free survival: 171 of 252 significant genes were found with univariate hazard ratio over 1 (Supplemental Table 1). With an independent text mining application querying NCBI database: GeneValorization [16], it could be able to establish a Ferroptosis/ECM remodeling molecular score in basal breast cancer [18] based on the expression of eleven related genes: TNF, IL6, SET, CDKN2A, EGFR, HMGB1, KRAS, MET, LCN2, HIF1A, TLR4. Except CDKN2A and TNF, the majority of these implicated genes have been validated as been expressed in TNBC cellular model MDA-MB-231 and regulated by ferroptosis activators: erastin and RSL3 [8,9].

HIF-1α is an important regulator of the lipid metabolism [29]. Hypoxia-induced lipid metabolism reprogramming results in fatty acid accumulation, which promotes tumor growth and survival upon reoxidation [30]. Ferroptosis induction with erastin or RSL3 showed a significant up regulation of HIF1A in TNBC cell line. In the context of increased hypoxia/HIF1A and ECM stiffness in chemoresistant tumors, high expression of HIF1A could be adverse because its is leading to upregulation of ITGA5, activation of the downstream FAK/Src signaling pathways and repression of miR-326 which targets fibronectin (FN1), an extracellular matrix (ECM) central chemoresistance driver gene [31]. The adipokine Lipocalin-2 (LCN2) has been demonstrated to be an ECM regulator through its association with the ECM protease matrix metalloproteinase-9 (MMP-9) [32]. It has been shown that LCN2 knockout in human breast cancer cell line MDA-MB-231 ameliorates erastin-mediated ferroptosis and increases cisplatin vulnerability [33].

Adipocytes constitute the main cell component of ECM in breast cancer [34]. Cancer-associated adipocytes (CAAs) are localized at the invasive front of breast tumor and exhibit a modified phenotype, a loss of lipid content, a decrease in late adipocyte differentiation markers and overexpression of inflammatory cytokines and proteases [35]. IL6 up regulation was promoted by ferroptosis activators in MDA-MB-231 TNBC cellular model (Figure 6D). In breast tumors, IL6 is secreted by CAAs which play essential roles in favor of proliferation, angiogenesis, dissemination, invasion and metastasis of breast cancer [36]. In breast cancer, the pleiotropic cytokine IL-6 production is related to the development of stem cell phenotype, angiogenesis, cachexia, and resistance to therapy [37].

EGFR overexpression was found under ferroptosis stimulation in MDA-MB-231 cells (Figure 6D). EGFR promoted TNBC cell clustering and blockade of EGFR successfully abolished tumor cell cluster formation [38]. It has been shown that inhibition of EGFR signaling pathway significantly suppressed cell viability of TNBC cells and reduced fraction of CSCs with intracellular enhancement of lipid peroxidation when TNBC cells exposed to erastin [39].

MET is known to be implicated in chemotherapy resistance especially resistance to targeted therapies, including those targeting EGFR, BRAF and MEK but also contributes to cytotoxic chemotherapy resistance [40]. Its ligand, HGF is a pleiotropic factor produced by mesenchymal cells in the stroma and as such is widely distributed in the extracellular matrix of most tissues [41]. Dysregulation of the MET/HGF pathway leads to uncontrolled cell proliferation and oncogenesis and is observed in multiple tumour types [42]. There is few link between MET and ferroptosis except that HGF exacerbated pancreatic cancer cell ferroptosis resistance [43]. In the present work, it was observed that RSL3 ferroptosis activator induced overexpression of MET in MDA-MB-231 TNBC cellular model.

In MDA-MB-231, it was observed a significant down regulation of High-mobility group box 1 (HMGB1). HMGB1 control chromosomal structure and is considered as a redox-sensitive protein. HMGB1 is implicated in regulating stress responses to oxidative damage and cell death, which are closely related to the pathology of inflammatory diseases such as cancer. HMGB1 can also be released into the extracellular space and function as a damage-associated molecular pattern protein during ferroptosis [44]. HMGB1 is known to act via the NRF2/GPX4 axis to repress ferroptosis in mesangial cells in response to high glucose [45].

In ER-positive breast cancer, SET overexpression reduced tamoxifen-induced antitumor effects and may contribute to the failure of the tamoxifen treatment by modulating the ER signaling [46]. In MDA-MB-231 TNBC cellular model, ferroptosis activation shown a major repression of SET proto-oncogene.

Ferroptosis activators induced repression of KRAS oncogene expression in MDA-MB-231 TNBC cellular model. KRAS mutations are known as very infrequent in triple-negative breast tumors [47], but in basal breast cancer, KRAS have been shown to promote the mesenchymal features of this aggressive cancer [48].

Overexpression of Toll-like receptor-4 (TLR4) in human tumors often correlates with chemoresistance and metastasis. Depletion of TLR4 in naturally overexpressing MDA-MB-231 cells downregulated prosurvival genes concomitant with 2- to 3-fold reduced IC(50) to paclitaxel in vitro and a 6-fold decrease in recurrence rate in vivo [49]. Effect of ferroptosis activators Erastin and RSL3 on MDA-MB-231 showed an opposite effect on TLR4 expression.

In the present work, by a text mining approach integrated in distinct transcriptome dataset, we could elucidate an adverse molecular program which crosstalk between ferroptosis and extracellular remodeling functionalities during basal breast cancer. This adverse regulated program allowed to compute a molecular expression score that could be promising to evaluate response of ferroptosis target therapies in breast cancer.