A systems-level analysis of the mutually antagonistic roles of RKIP and BACH1 in dynamics of cancer cell plasticity

A systems-level analysis of the mutually antagonistic roles of RKIP and BACH1 in dynamics of cancer cell plasticity | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results A systems-level analysis of the mutually antagonistic roles of RKIP and BACH1 in dynamics of cancer cell plasticity Sai Shyam , R Soundharya , Manas Sehgal , View ORCID Profile Mohit Kumar Jolly doi: https://doi.org/10.1101/2023.07.08.547338 Sai Shyam 1 Centre for BioSystems Science and Engineering, Indian Institute of Science , Bangalore, 560012, India Find this author on Google Scholar Find this author on PubMed Search for this author on this site R Soundharya 1 Centre for BioSystems Science and Engineering, Indian Institute of Science , Bangalore, 560012, India Find this author on Google Scholar Find this author on PubMed Search for this author on this site Manas Sehgal 1 Centre for BioSystems Science and Engineering, Indian Institute of Science , Bangalore, 560012, India Find this author on Google Scholar Find this author on PubMed Search for this author on this site Mohit Kumar Jolly 1 Centre for BioSystems Science and Engineering, Indian Institute of Science , Bangalore, 560012, India Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Mohit Kumar Jolly For correspondence: mkjolly{at}iisc.ac.in Abstract Full Text Info/History Metrics Supplementary material Preview PDF Abstract Phenotypic plasticity is a hallmark of cancer metastasis. Epithelial-mesenchymal transition (EMT) is an important axis of phenotypic plasticity. Raf kinase-B inhibitor protein (RKIP) and BTB and CNC homology 1 (BACH1) are two proteins reported to influence EMT. In breast cancer, they act antagonistically, but the exact nature of their roles in mediating EMT and associated other axes of plasticity remains unclear. Here, analysing transcriptomic data, we reveal their antagonistic trends in a pan-cancer manner, in terms of association with EMT, metabolic reprogramming and immune evasion via PD-L1. Next, we developed and simulated a mechanism-based gene regulatory network that captures how RKIP and BACH1 engage in feedback loops with drivers of EMT and stemness. We found that RKIP and BACH1 belong to two separate “teams” of players – while BACH1 belonged to the one that drove pro-EMT, stem-like and therapy-resistant cell-states, RKIP is a member of a team that enables pro-epithelial, less stem-like and therapy-sensitive phenotypes. Finally, we observed that low RKIP levels and concomitant upregulated BACH1 levels associated with worse clinical outcomes in many cancer types. Together, our systems-level analysis indicates that the emergent dynamics of underlying regulatory network underlie the antagonistic patterns of RKIP and BACH1 with various axes of cancer cell plasticity, as well as with patient survival data. Introduction The process of cancer metastasis is driven by phenotypic plasticity, i.e. dynamic and reversible adaptation of disseminating cancer cells to different microenvironments that they encounter along the journey. Phenotypic switching among the epithelial (E), mesenchymal (M) and hybrid E/M phenotypes through Epithelial to Mesenchymal Transition (EMT) and its reverse Mesenchymal to Epithelial Transition (MET) constitute an important axis of phenotypic plasticity during metastasis ( Brabletz et al., 2018 ; Jia et al., 2019 ). Metabolic reprogramming is another key axis of plasticity and a hallmark of cancer metastasis ( Celià-Terrassa and Kang, 2016 ). Different axes of plasticity – EMT/MET, metabolic switching and immune evasion – are often interconnected, thus enabling cancer cells to exist in distinct cell-states ( Dongre et al., 2017 ; Jia et al., 2021 ; Noman et al., 2017 ), and promoting phenotypic plasticity and consequent non-genetic heterogeneity. Understanding the dynamics of the interconnected axes of plasticity is critical to restrict metastasis. For instance, EMT often leads to increased levels of PD-L1 – a transmembrane molecule leading to immune escape of cancer cells (Chen et al., 2014). Further, knockdown of PD-L1 could reverse EMT ( Alsuliman et al., 2015 ). Similarly, EMT and tamoxifen resistance in ER+ breast cancer cells can drive each other ( Sahoo et al., 2021 ). Such bidirectional connections are often mediated by multiple feedback loops among the molecules driving cell plasticity along these multiple axes. Breaking the miR-200/ZEB1 mutually inhibitory feedback loop in breast cancer cells through CRISPR/Cas9 can reduce cancer cell metastasis ( Celià-Terrassa et al., 2018 ). Thus, mapping the different feedback loops that can govern cell plasticity is of fundamental importance. Recent reports have identified a mutually inhibitory feedback loop between the Raf kinase inhibitor protein (RKIP) and BTB and CNC homology 1 (BACH1) in regulation of EMT and metastasis in breast cancer (Lee et al., 2014; Wan et al., 2023 ). BACH1 represses RKIP transcriptionally and RKIP can inhibit BACH1 via microRNA let-7. RKIP is a metastasis suppressor that acts along the RAF1/MEK/ERK pathway to regulate cell proliferation and migration ( Dong et al., 2021 ). Its exogenous expression in metastatic breast cancer cells can suppress invasion, intravasation and metastasis in xenograft mouse models ( Dangi-Garimella et al., 2009 ). On the other hand, BACH1 can promote breast cancer metastasis ( Yun et al., 2011 ). While their antagonistic roles in breast cancer have been reported, it remains unclear whether this antagonism is seen in other cancers, and how do these molecules regulate different axes of cellular plasticity implicated in metastasis. Here, we first investigated whether RKIP and BACH1 show antagonistic trends across different cancer types using transcriptomic data from The Cancer Genome Atlas (TCGA). We found RKIP and BACH1 to be anti-correlated with each other in majority of the cancer types. Moreover, while BACH1 correlated positively with EMT and PD-L1 but negatively with oxidative phosphorylation and fatty acid oxidation, RKIP showed opposite trends in a pan-cancer manner. These trends were recapitulated in ER+ breast cancer datasets, where BACH1 also correlated negatively with ESR1 (Estrogen Receptor), but RKIP correlated positively with it. To better understand these consistent patterns in transcriptomic signatures, we constructed and simulated a mechanism-based gene regulatory network (GRN) that incorporated the feedback loops formed among RKIP, BACH1 and other master regulators of cancer cell plasticity such as ZEB1, miR-200, LIN28 and let-7. Our analysis of GRN identified that RKIP and BACH1 belonged to two mutually repressing “teams” of players – one that was comprised of pro-EMT (ZEB1, LIN28, SNAIL, SLUG) players, and the other constituted pro-MET (miR-200, let-7, miR-145, CDH1) ones. The existence of these “teams” enables the BACH1-high cells to display a hybrid E/M and mesenchymal state exhibiting a stem-like behaviour, as well as opposite trends in terms of association with patient survival across cancer types. Together, our results explain the emergent dynamics of underlying GRN that can underlie the observed antagonistic behaviour of RKIP and BACH1 in a pan-cancer manner. Results RKIP and BACH1 display mutual antagonistic patterns across many cancer types in TCGA RKIP and BACH1 have been reported as mutually inhibitory players in breast cancer; while RKIP is anti-metastatic, BACH1 is pro-metastatic (Lee et al., 2014). To investigate whether they show antagonistic trends consistently across other cancers, we analysed TCGA data from 35 different cancer types. RKIP and BACH1 were found to be negatively correlated with each other in 31 out of 35 (88.57%) cancer types ( Fig 1A ). To understand their association with EMT, we calculated the ssGSEA score for each sample for the KS-Epithelial (KS-Epi) and KS-Mesenchymal (KS-Mes) signatures. KS-Epi signature comprises of genes that are upregulated in epithelial cells in a pan-cancer manner, while genes in KS-Mes signatures are upregulated in mesenchymal cells ( Tan et al., 2014 ). The ssGSEA score of KS-Epi signature correlated positively with RKIP levels in only 15 cancer types and negatively with BACH1 expression in 14 cancer types, indicating a rather ambivalent association. On the other hand, the KS-Mes scores were negatively correlated with RKIP levels in 57.14% (20 out of 35) cancer types, but positively correlated with BACH1 expression in 88.57% (31 out of 35) of cases ( Fig 1A ). Further, RKIP showed a negative correlation with KS EMT score in 57.14% (20 out of 35) cancers but BACH1 showed positive correlation with it in 77.14% (27 out of 35) cases. The higher the KS EMT score, the more mesenchymal the sample is ( Chakraborty et al., 2020 ). Together, these results suggest that RKIP and BACH1 show antagonistic trends in a pan-cancer manner, with RKIP associating with an epithelial phenotype while BACH1 with a mesenchymal phenotype. Download figure Open in new tab Figure 1: Correlation of RKIP and BACH1 expression with cellular plasticity signatures in TCGA. Heatmaps representing Spearman’s correlation coefficient of RKIP and BACH with A) each other, KS Score and ssGSEA scores of KS-Epi and KS-Mes gene signatures, B) ssGSEA scores of Hallmark gene signatures (GS) for OXPHOS, Glycolysis and FAO, C) ssGSEA scores for PDL1 gene signature (GS) and D) ssGSEA scores for ferroptosis gene signature (GS). * denotes p-value < 0.05. During EMT, cells can also exhibit metabolic plasticity, typically leading to decreased oxidative phosphorylation and fatty acid oxidation but increased glycolysis ( Muralidharan et al., 2022a ). To investigate how RKIP and BACH1 expression correlates with these pathways, we calculated the ssGSEA scores for OXPHOS, FAO and glycolysis gene signatures in TCGA samples. OXPHOS ssGSEA scores correlated positively with RKIP in 88.57% (31 out of 35) of the cancer types while BACH1 correlated negatively in 91.14% (32 out of 35) of the cancer types ( Fig 1B ). Similar trends were observed with FAO, with RKIP correlating positively and BACH1 correlating negatively in 31 out of 35 cancer types. On the other hand, ssGSEA scores for glycolysis signature did not show such strong trends – they correlated negatively with RKIP in 14 cancer types, and positively with BACH1 in 9 of them. Also, while FAO and OXPHOS ssGSEA scores correlated positively with each other, glycolysis did not show consistent negative correlations with either of them ( Fig S1A ) . Overall, RKIP and BACH1 exhibited opposite patterns in terms of their association with OXPHOS and FAO. Another axis of phenotypic plasticity coupled with EMT is immune evasion. As cells undergo a partial or complete EMT, they display increased levels of PD-L1, an immune checkpoint molecule that evades attack by the immune system ( Dongre et al., 2021 , 2017 ). Thus, we probed how the expression levels of PD-L1 and scores for PD-L1 associated gene signature correlate with RKIP and BACH1. We observed that among 35 cancer types, PD-L1 expression correlated negatively with RKIP expression in 22 of them, and positively with BACH1 in 32 of them ( Fig S1B ). Similarly, the ssGSEA scores of PD-L1 gene signature correlated negatively with RKIP in 68.57% (24 out of 35) but positively with BACH1 in 88.57% (31 out of 35) of the cancer types ( Fig 1C ). Further, cells undergoing EMT have been shown to be vulnerable to ferroptosis, an iron-dependent cell-death program ( Viswanathan et al., 2017 ). To examine how RKIP and BACH1 associate with ferroptosis, we calculate the ssGSEA scores of a ferroptosis-based gene signature in cancer ( Lu et al., 2021 ). Consistently, we found that RKIP associates negatively with it, but BACH1 ( Fig 1D ). Overall, BACH1 is likely associated with a more mesenchymal, glycolytic, ferroptosis-sensitive and immune-evasive phenotype, but RKIP tends to promote an epithelial, immune-sensitive and ferroptosis-insensitive cell-state dependent on OXPHOS and FAO. Association of RKIP and BACH1 with EMT and tamoxifen resistance in ER+ breast cancer Next, we focused on breast cancer, given the earlier observations about mutually antagonistic roles of RKIP and BACH1 in breast cancer ( Lee et al., 2019 , 2014 ; Yesilkanal et al., 2021 ). We first analysed how RKIP and BACH1 correlate with ESR1 – gene encoding for Estrogen Receptor alpha (ERα). Higher levels of ERα are often associated with improved response to anti-estrogen therapies such as tamoxifen, and better patient survival ( Burns and Korach, 2012 ). In TCGA breast cancer samples, we observed that ESR1 correlated positively with RKIP (ρ = 0.11), but negatively with BACH1 (ρ = − 0.29) ( Fig 2A ). ZEB1, an EMT-inducing transcription factor (EMT-TF), has been shown to hypermethylated the ERα promoter and confer tamoxifen resistance ( Zhang et al., 2017 ). Thus, we evaluated association of RKIP and BACH1 with ZEB1. ZEB1 correlated significantly positively with BACH1 (ρ = 0.52) and negatively with RKIP (ρ = − 0.32) ( Fig 2B ). Similarly, BACH1 correlated positively with another EMT-TF (SNAI2) and an EMT marker vimentin (VIM) but negatively with OVOL2, an MET-inducing transcription factor ( Saxena et al., 2022 ), showing opposite trends as those seen for RKIP ( Fig S2A ). These results indicate that BACH1 associates with EMT, but RKIP associates with MET in breast cancer. To investigate the associations of RKIP and BACH1 with metastasis, we used their respective pathway metastasis signatures ( Lee et al., 2013 ; Yun et al., 2011 ). In TCGA breast cancer samples ( Fig 2C ) and in other TCGA cancer types ( Fig S2B ), ssGSEA scores of RKIP Pathway Metastasis Signature (RPMS) correlated negatively with RKIP, while those of BACH1 Pathway Metastasis Signature (BPMS) correlated positively with BACH1. Consistently, RPMS and BPMS ssGSEA scores correlated positively with one another in a pan-cancer manner ( Fig S2B ), thus endorsing the pro-metastatic role of BACH1 and anti-metastatic role of RKIP. Download figure Open in new tab Figure 2: Antagonistic trends of RKIP and BACH1 with EMT & associated axes in breast cancer. A, B) Scatter plots showing correlations of RKIP and BACH1 with ESR1, ZEB1 in TCGA breast cancer samples. C, D) Same as A, B but for RKIP pathway metastasis signature (RPMS), BACH1 pathway metastasis signature (BPMS). E) Heatmap depicting Spearman’s correlation coefficient in 6 ER +ve Breast cancer datasets - GSE17705, GSE24202, GSE43495, GSE6532, GSE67916, GSE9195. i) Correlation of RKIP and BACH1 with BPMS, RPMS; ii) Correlation with ferroptosis and PD-L1 gene signatures (GS). iii) correlation of BPMS and RPMS with metabolic axes - OXPHOS, FAO, glycolysis. In these heatmaps, ** : p-value < 0.05, * : p-value< 0.1 (Spearman’s correlation analysis). We next focused on six ER+ breast cancer datasets that we had previously analysed from the perspective of EMT and tamoxifen resistance driving each other in ER+ breast cancer ( Sahoo et al., 2021 ): GSE6532, GSE9195, GSE17705, GSE24202, GSE43495 and GSE67916. We noticed resonating trends in these datasets as noted in TCGA samples: RKIP inversely correlated with ssGSEA scores of RPMS and ferroptosis gene signature, while BACH1 correlated positively with ssGSEA scores of BPMS and ferroptosis signature ( Fig 2E, i-ii ). Further, both RPMS and BPMS scores both correlate positively with those of PD-L1 signature but negatively with OXPHOS ( Fig 2E, ii-iii ), reminiscent of observations of PD-L1 activity being anti-correlated with OXPHOS across carcinomas ( Muralidharan et al., 2022b ). Across these datasets, RPMS correlates positively with glycolysis but negatively with FAO, while BPMS shows relatively weaker trends ( Fig 2E, iii ). This analysis further strengthens the pan-cancer trends that we noticed earlier in terms of antagonistic association of RKIP and BACH1 with multiple axes of phenotypic plasticity. Underlying gene regulatory networks reveal RKIP and BACH1 as members of mutually antagonistic “teams” To understand the mechanisms that can explain the observations about RKIP and BACH1 showing opposite trends with respect to metastatic propensity, we identified a minimal core underlying gene regulatory network (GRN) in breast cancer that incorporates the feedback loops that RKIP and BACH1 are involved in with key players of EMT, tamoxifen resistance and stemness, building on our previous efforts to connect these axes of plasticity ( Pasani et al., 2021 ; Sahoo et al., 2021 ). This network is not meant to be exhaustive in terms of connections RKIP and BACH1 have with players controlling these axes of plasticity, but demonstrate a core network structure that may be sufficient to explain the correlation-based observations of RKIP and BACH1 noted earlier. Broadly speaking, this network has three core modules: EMT, tamoxifen resistance and stemness ( Fig 3A ). The EMT module comprises EMT-TFs SNAIL, SLUG and ZEB1, EMT-inhibiting micro-RNA-200 family and E-cadherin ( CDH1 ), a key cell-cell adhesion molecule that maintains tight junctions among epithelial cells ( Sahoo et al., 2021 ). The stemness module is composed of OCT4, LIN28, let-7 and miR-145. LIN28 and let-7 engage in a mutually inhibitory loop, and so do miR-145 and OCT4 ( Jolly et al., 2016 ; Pasani et al., 2021 ). In tamoxifen resistance module, we include ERα66 and ERα36, two variants of estrogen receptor ( ESR1 ) – ERα66 associates with tamoxifen-sensitive cell-state, while elevated levels of ERα36 drive resistance to tamoxifen in breast cancer cells ( Shi et al., 2009 ). ERα66 can suppress ERα36 levels ( Zou et al., 2009 ) directly, while ERα36 can inhibit ERα66 by activating ZEB1 ( Sahoo et al., 2021 ). Other links across modules involve inhibition of LIN28 by miR-200 ( Kong et al., 2010 ) and inhibition of ERα66 by miR-145 ( Spizzo et al., 2009 ) as well as mutual inhibition between miR-145 and ZEB1 ( Jolly et al., 2016 ), and that between SLUG and ERα66 (Sahoo et al., 2021a). RKIP and BACH1 associate with these modules through the following links: a) BACH1 can self-inhibit and activate SLUG ( Igarashi et al., 2021 ), b) RKIP and SNAIL inhibit each other ( Beach et al., 2007 ; Shvartsur et al., 2017 ), and c) BACH1 represses RKIP directly, while RKIP inhibits BACH1 via let-7 (Lee et al., 2014). Download figure Open in new tab Figure 3: Gene regulatory network underlying antagonistic roles of RKIP and BACH1 as members of opposite “teams”. A) Regulatory network showing interconnections of RKIP and BACH1 with the modules of EMT, tamoxifen resistance and stemness. B) Adjacency Matrix for the network shown in A. C) Influence Matrix for path length = 10 for network shown in A. D) Network randomisation schematic (inset) & histogram of team strength of 1000 random networks generated. The red line indicates the team strength of the wild-type (WT) network, i.e. shown in A. E) Scatterplots of Principal Component 1 (capturing 48.8% variance) and Principal Component 2 (capturing 19.5% variance) shaded by BACH1 and RKIP Levels. F) Bar plot showing the PC1 loading coefficients of different genes in the network. The bar colors show the same antagonistic team behavior of the teams identified in C. Before simulating the emergent dynamics of this network, we generated corresponding adjacency matrix. This matrix – showing only the direct links (both activation and inhibition) between different nodes in a network – is rather spare ( Fig 3B ). Because the influence of one node on another can also be mediated through indirect links, we derived the influence matrix for this network ( Hari et al., 2022 ), up to a path length of eight edges, and performed hierarchical clustering ( Fig 3C ). The clustering revealed the presence of two “teams” such that members within a team effectively activated one another, while members across two teams effectively inhibited each other. One team consisted of players corresponding to an epithelial phenotype (let-7, miR-200, CDH1) while the other team comprised drivers of EMT (ZEB, SNAIL, SLUG) and stemness (LIN28, OCT4). As expected, ERα66 was found to be a part of the MET-promoting “team”, while ERα36 belonged to an EMT-driving one. Similarly, RKIP and BACH1 were found in two opposite “teams”, consistent with our previous TCGA analysis showing RKIP to be pro-epithelial but RKIP to be pro-EMT. Next, we quantified the team strength for this influence matrix and found it to be 0.37 (on a scale of 0 to 1). To examine whether this observation of “teams” is unique to this network, we generated 1000 random networks by shuffling the edges among the nodes (following the randomization schematic; Fig 3D inset) and determined their corresponding team strength. We observed the team strength of the wild type (WT) network to be much greater than that of all 1000 randomly generated networks ( Fig 3D ). This observation suggests that the organization of these molecular players into “teams” in this network is not a co-incidence, rather a unique topological signature that can possibly facilitate specific couplings between cellular behaviour such as association of EMT with stemness and tamoxifen resistance ( Jolly and Celia-Terrassa, 2019 ; Wang et al., 2019 ) and opposite roles of RKIP and BACH1 in EMT ( Cessna et al., 2022 ; Igarashi et al., 2021 ). Further, we simulated the dynamics of this GRN by representing the regulatory interactions through a set of coupled ordinary differential equations (ODEs) across an ensemble of kinetic parameters and initial conditions, using a tool called RACIPE ( Huang et al., 2018 ). The output of RACIPE is the set of steady-state solutions obtained which are then z-normalized for a better comparison across the expression patterns. Principal component analysis (PCA) performed on RACIPE solutions indicated that over 60% of the variance could be captured by first 2 principal components (PCs) ( Fig S3A ) - PC1 accounted for 48.8%, while PC2 accounted for 19.5%. Colouring the PCA plots based on z-normalized levels of RKIP and BACH1 revealed that clusters showing higher BACH1 levels had lower RKIP levels and vice versa ( Fig 3E ). Next, we plotted the loading coefficients of each node in the network along the PC1 axis to understand their contributions to PC1 variance. We observed that genes earlier identified to be a part of the epithelial team from the influence matrix analysis had PC1 coefficients less than zero, while all of them in mesenchymal team (with SNAIL being the only exception) had these coefficients greater than zero. This analysis reveals that PC1 is largely able to recapitulate the members belonging to two different “teams” identified via influence matrix. The existence and constitution of these “teams” is further validated by largely bimodal distribution of all network nodes, PCA correlation circle and pairwise correlation matrices showcasing positive correlation among members within a team, and negative across “teams” ( Fig S3B-C ). Regulatory network dynamics explains the association of BACH1 with stem-like cell-state Next, we investigated whether this regulatory network is capable of multi-stable behaviour, i.e. allowing for the co-existence of multiple cell-states that can reversibly switch among themselves. We segregated the parameter sets generated by RACIPE based on their corresponding number of stable states ( Fig 4A ). Our analysis revealed that only 4.58% of solutions correspond to a mono-stable state; the remaining 95.42% solutions associated with two or more co-existing cell-states, highlighting the underlying multistable dynamics of this GRN that can support phenotypic plasticity. Download figure Open in new tab Figure 4: Association of RKIP and BACH1 with EMT and stemness states. A) Pie chart depicting the number of parameter sets giving rise to different number of stable states. Data from 3 independent replicate RACIPE simulations shown as mean ± standard deviation). B) Heatmap of stable steady-state solutions for network shown in Fig 3A , simulated from RACIPE. Color bar represents the relative levels of individual genes (z-normalized and log2 expression values). C) i) Kernel Density Plot of EM Score (= ZEB +SLUG – miR200 – CDH1). Red dotted lines indicate the segregation of Epithelial, Hybrid and Mesenchymal phenotypes based on minima of the fitted Gaussians. ii) Kernel Density Plot of SN Score (= OCT4 +LIN28 – mir145 – let7). Red dotted lines indicate the segregation of stem-like and non-stem-like phenotypes based on minima of the fitted Gaussians. D) Scatterplot of all RACIPE solutions, on SN Score vs EM Score axes. Segregation done based on thresholds identified in C i, ii. E) i) Table depicting the conditional probability percentage of gene expression status given the epithelial, mesenchymal or hybrid Phenotypes. ii) Same as i) but given the stem-like and non-stem-like categorization. BACH1 (or RKIP) +ve indicates z-score expression of BACH1 (or RKIP) > 0, and BACH1 –ve indicates it < 0. Further, we plotted the ensemble of stable states obtained from RACIPE as a heatmap ( Fig 4B ). It showed expected clustering patterns and the dominance of two distinct states – (high LIN28, high OCT4, high BACH1, high ERα36, high SLUG, high ZEB, low let-7, low ERα66, low SNAIL, low miR-145, low miR-200, low CDH1) and (low LIN28, low OCT4, low BACH1, low ERα36, low SLUG, low ZEB, high let-7, high ERα66, high SNAIL, high miR-145, high miR-200, high CDH1), thus recapitulating the “teams” seen in influence matrix of the GRN. These two states correspond to a mesenchymal, stem-like, tamoxifen-resistant and epithelial, non-stem-like, tamoxifen-sensitive state respectively. Interestingly, higher levels of RKIP were seen in most RACIPE solutions that correspond to an epithelial cell-state, as well as a subpopulation of mesenchymal cell-state subpopulation. This observation is consistent with a lower magnitude of loading coefficient of RKIP in explaining PC1 variance as compared to that of BACH1 ( Fig 3D ). Thus, our simulations suggest that BACH1 associates strongly with a partial/full EMT state, relative to the association of RKIP with an epithelial one. To better elucidate the functional mapping of phenotypes from the RACIPE expression data, we defined an Epithelial Mesenchymal (EM) score as EM Score = (ZEB1 + SLUG − miR200 − CDH1)/4. ZEB1 and SLUG are characteristic EMT markers, while miR-200 and CDH1 are those for an epithelial state. The higher the EM score, the more mesenchymal the corresponding cell-state is. A histogram of EM scores thus obtained across all RACIPE solutions indicated three distinct populations which can be characterized as epithelial, mesenchymal and hybrid ones ( Fig 4C,i ). Similarly, to quantify plasticity along the stemness axes, we defined a stemness (SN) score as SN score = (OCT4 + LIN28 − miR145 − let7)/4. Extremely high or low levels of OCT4 and LIN28 are associated with non-stem-like states ( Jolly et al., 2014 ; Karwacki-Neisius et al., 2013 ; Niwa et al., 2000 ), thus we ascribe intermediate levels of SN scores to a stem-like state, while very high or low levels of it to a non-stem-like state, based on the histogram ( Fig 4C,ii ). The z-normalized RACIPE solutions were then projected on a scatter plot with corresponding SN and EM scores as the axes. While the scatter plot showed a general positive correlation between EM and SN scores (ρ = 0.701), we noticed associations of all the three phenotypes along the EMT axis with both stem-like and non-stem-like states ( Fig 4D ). To better understand how RKIP and BACH1 levels associated with these axes, we divided the RACIPE steady state solution ensemble into BACH1 (+ve) and BACH1 (-ve) (and similarly for RKIP (+ve) and RKIP (-ve)) cases based on their corresponding z-normalized values. We observed that 63.48% of in silico cells (considering each steady-state solution of RACIPE as equivalent to a cell in an experimental population-level setting) showing a hybrid E/M phenotype were BACH1(+ve) ( Fig 4E,i ). Similarly, a mesenchymal cell was three times more likely to be BACH1(+ve) as compared to being BACH1(-ve) (74.94% vs 25.06% cases) ( Fig 4E,i ). Consistent with experimental observations of partial and/or full EMT phenotypes with stemness ( Brown et al., 2022 ; Lourenco et al., 2020 ; Lüönd et al., 2021 ), we found a stem-like cell to be 2.3 times (69.64% vs 30.36%) more likely to be BACH1(+ve) as compared to being BACH1 (-ve) ( Fig 4E, ii ). On the other hand, the RKIP(+ve) or RKIP(-ve) cells were not enriched in any subpopulation along the EMT or stemness axes. Together, we can infer that BACH1 associates with a hybrid E/M and/or mesenchymal stem-like phenotype. High BACH1 expression levels associate with worse patient survival in many cancer types To identify how RKIP and BACH1 expression levels impact survival outcome, we computed the hazard ratios for different combinations of RKIP and BACH1 (RKIP+ BACH1-, RKIP+ BACH1+, RKIP-BACH1-) with RKIP-BACH1+ as reference. High BACH1 expression along with low RKIP expression led to significantly worse overall survival when compared to high RKIP and low BACH1 samples in three different cancers: lung adenocarcinoma (LUAD), liver hepatocellular carcinoma (LIHC) and pancreatic adenocarcinoma (PAAD). Similar trends were recapitulated for disease specific survival, disease-free survival and progression-free survival in these cancers. Conversely, RKIP+ BACH1+ samples had significantly better progression free survival compared to the RKIP-BACH1+ in LIHC, indicating that higher RKIP expression may be linked to better survival outcomes ( Fig 5 ). Consistent trends are also observed in mesothelioma, breast cancer, renal clear cell carcinoma, uveal melanoma and thyroid carcinoma ( Fig S4 ). These pan-cancer observations carried out consistently indicate that high BACH1 levels are linked to worse survival. Download figure Open in new tab Figure 5: Forest plots showing antagonistic trends of RKIP and BACH1 with varying survival types across different cancers. Reference set of TCGA sample considered here is R-B+ (RKIP-low, BACH1-high). p-values are based on log-rank test, and those with significant differences p < 0.05, p<0.01 and p<0.001 are marked with *, **, and ***, respectively. A) Hazard Ratios (HR) for Overall Survival in PAAD (left), LUAD (middle) and LIHC (right). B) HR for Disease Specific Survival in PAAD (left), LUAD (middle) and LIHC (right). C) HR for Disease-Free Survival in in PAAD (left), LUAD (middle) and LIHC (right). D) HR for Progression -Free Survival in PAAD (left), LUAD (middle) and LIHC (right). Discussion Cellular decision-making is often governed by mutually inhibitory feedback loops, also called as a ‘toggle switch’. Such network motifs are commonly reported in various decision-making contexts, such as GATA1-PU.1 toggle switch determining differentiation of a common myeloid progenitor into erythroid cells (GATA1 >> PU.1) or myeloid cells (PU.1 >> GATA1) ( Zhou and Huang, 2011 ). Similar antagonism has been witnessed among many EMT and MET driving transcription factors (TFs) such that the MET-TFs (GRHL2, KLF4, OVOL2 etc.) and EMT-TFs (ZEB1, SNAIL, TWIST etc.) mutually repress one another ( Chung et al., 2016 ; Roca et al., 2013 ; Subbalakshmi et al., 2021 ). Such mutually antagonistic “teams” of players – also seen in small cell lung cancer GRN ( Chauhan et al., 2021 ) – can thus contribute to cellular plasticity ( Hari et al., 2022 ). Here, we show that RKIP and BACH1 belong to two antagonistic “teams” – while RKIP belongs to a “team” promoting an epithelial and tamoxifen-sensitive cell-state in breast cancer, BACH1 is a part of “team” enabling a mesenchymal and tamoxifen-resistant phenotype. While the antagonistic role of RKIP and BACH1 in EMT has been well-reported, their potential involvement in tamoxifen resistance is only beginning to be investigated through mediators such as TANK-binding kinase 1 (TBK1) ( Liu et al., 2022 ; Wei et al., 2014 ). Consistently, their contrasting roles are emerging in the contexts of ferroptosis ( Wenzel et al., 2017 ; Xie et al., 2023 ) and metabolic reprogramming ( Lee et al., 2019 ) in specific cancer types. Particularly, activation of NRF2 – a master regulator of anti-oxidant program in cells – stabilizes BACH1 ( Lignitto et al., 2019 ). NRF2 activity can be enhanced by loss of RKIP ( Al-Mulla et al., 2012 ), enabling another ‘toggle switch’ between RKIP and BACH1. NRF2 has also been reported to maintain cells in a hybrid epithelial/mesenchymal (E/M) phenotype ( Vilchez Mercedes et al., 2022 ) and enhance stemness and chemoresistance ( Al-Mulla et al., 2012 ; Kim et al., 2018 ). Thus, our results about association of BACH1 with a hybrid E/M stem-like state unify the previous observations about role of BACH1 in controlling multiple axes of plasticity. The contrasting roles of RKIP and BACH1 in mediating stemness/dedifferentiation lends further credence to our model simulations. BACH1 can activate CD44 and MAPK signaling in lung cancer stem cells (CSCs) and stimulate lung cancer metastasis; its loss represses metastasis in xenograft models ( Jiang et al., 2021 ; Wiel et al., 2019 ). Similarly, BACH1 represses mesendodermal differentiation in embryonic stem cells, maintaining stem-cell identity ( Wei et al., 2019 ). Conversely, RKIP can repress NANOG in primary melanocytes, maintaining their differentiation state ( Penas et al., 2020 ). Chemical induction of RKIP can degrade SOX2, inhibit tumor growth and promote differentiation of schwannoma into mature Schwann cells ( Cho et al., 2022 ). These observations, together with the contrasting roles of RKIP and BACH1 on tumor cell migration ( Davudian et al., 2016 ; Zhu et al., 2018 ), can explain the association of higher BACH1 levels with enhanced metastasis and poor patient prognosis ( Chen et al., 2023 ; Han et al., 2019 ; Wiel et al., 2019 ). Our results demonstrate that these trends of association of RKIP and BACH1 with clinical outcomes are largely consistent across cancer types as well as survival metrics (overall survival, progression-free survival, disease-free survival), highlighting their behavior as members of two “teams” playing a tug-of-war – a pro-metastatic one (EMT, stem-like, drug-resistant) and an anti-metastatic one. Overall, our pan-cancer systems-level analysis reveals that RKIP and BACH1 can control multiple axes of plasticity (EMT, metabolic reprogramming, stemness) together in opposite directions, thus explaining their association with patient survival seen across cancer types. Our mechanistic model highlights that such “teams” of players can be an important network motif in co-ordinating changes along more than one axes of plasticity together, such as EMT, stemness and metabolic switching. Our results suggest breaking such “teams” as a possible therapeutic avenue to reduce the fitness of metastasizing cells, by limiting their phenotypic plasticity trajectories. Conflict of Interest The authors declare no conflict of interest. Author contributions SS, SR and MS performed research and analysed data. SS prepared initial draft of the article, and all authors edited it. MKJ designed and supervised research and acquired funding. Code and Data Availability https://github.com/saishyam1/RKIP_BACH1_Data_Codes.git Materials & Methods 1. Transcriptomic datasets We downloaded microarray data from NCBI GEO (GSE17705, GSE24202, GSE43495, GSE6532, GSE67916, GSE9195) using the GEOquery package on R, and mapped the probes to respective genes using their corresponding annotation files. We normalized gene-wise expression matrices on log2 base, before further analysis. We downloaded TCGA gene expression data of 35 different cancer types, obtained from the UCSCXena browser. These pre-processed datasets were available in transcripts per million (TPM) format and were directly used for analysis. Survival data was analysed using data from TCGA. Relevant cancer samples were split into two groups based on median: high PEBP1 (gene name for RKIP) vs low PEBP1, and high BACH1 vs low BACH1. ‘Kaplan -Meier curves for overall survival were plotted using the plotter on the ‘KMPlotter’ website ( Nagy et al., 2021 ). Additionally, ‘coxph’ function in R package ‘survival’ was employed to determine the hazard ratio (HR) and confidence interval (95% CI) for TCGA cohorts, and heatmaps were made using ‘ggplot2’. Normalized Single Sample Gene Set Enrichment Analysis (ssGSEA) scores based on specific input gene signatures ( Table S1 ) were calculated for each sample using ‘GSEApy’ python package ( Subramanian et al., 2005 ). For correlation analysis between any two variables, Spearman’s correlation coefficient has been used, using ‘spearmanr’ function from ‘SciPy’ Python library. 2. EMT scoring metric KS score is a metric to quantify the extent of EMT based on expression levels of specific epithelial and mesenchymal markers ( Tan et al., 2014 ). It uses 2 gene signatures: KS-Epi (genes associated with an epithelial phenotype), KS-Mes (genes associated with a mesenchymal phenotype). It plots two cumulative distribution functions (CDFs) based on expression levels of genes in KS-Epi and KS-Mes signatures. Distance between the 2 CDFs is calculated for each value, and the maximum value is taken as the statistic. The KS score values lie in the range [-1, 1]. 3. Identification of “teams” of players An adjacency matrix (Adj) represents the topology of a GRN in the form of a matrix. Rows represent the source node for a particular link, while the columns represent the target nodes. An inhibitory link is represented (in blue) with value -1; an activation link is represented (in red) with value 1. Path-length is defined as the number of consecutive links in a path that connects a source node to its corresponding target. In a GRN, the nodes not only influence their direct target (path length = 1) but also other nodes indirectly (path length > 1). To take into consideration these interactions, we define an influence matrix that captures both the direct and indirect interactions between nodes in the network up to a defined path length ( Hari et al., 2022 ). where Adj l represents the adjacency matrix (Adj) multiplied with itself l times. Adj max represents the Adj matrix with all inhibition links replaced by activation links. Adj l / Adj l max represents the element-wise division of values in Adj l by values in Adj l max . A positive value indicates activation, and a negative influence indicates inhibition. The higher the value in the influence matrix, the higher the influence of that specific source gene on the target. Hierarchical clustering is performed on the influence matrix to identify the clusters of genes functioning similarly. The team Strength of each cluster T KL is given as Team Strength of the entire network T S is given as 4. RACIPE analysis Random Circuit Perturbation (RACIPE) is a tool used to simulate the dynamics of GRNs. It takes as input the topology of a GRN and generates an ensemble of kinetic models for the given GRN. Here, we simulated the network shown in Fig 3A – the input file to RACIPE is given as Table S2 . For each kinetic model generated for the input topology file, RACIPE samples many initial parameters from the designated range for each parameter. The expression levels of a node in a GRN is determined by a set of ordinary differential equations given below. Here, X i is concentration of gene product of gene node X part of GRN and is the rate of change of the gene expression with respect to time. is basal production rate; is basal degradation rate of gene product X i . H s represents shifted Hills function that models the activation and inhibition links towards this (X i ) gene node in the GRN. n ji is Hill’s function coefficient and is threshold value of the Hill equation. X j is concentration of gene product X j , where X j is a node either activating or inhibiting X i . Since our GRN had 13 nodes, i and j can take integer values between 1 and 13. λ ji is fold change parameter. Throughout this study, we simulated GRNs for 10,000 parameter sets. RACIPE gives log 2 normalized steady state gene expressions of each gene product as output. This output steady state gene expression data was then z-normalized and used for analysis. To check for bimodality in solutions obtained from RACIPE, we used Sarle’s bimodality coefficient: Here, ‘s’ is the skew of the distribution and ‘k’ is the kurtosis. BC can take values between 0 and 1, with BC values greater than 0.55 indicating bimodal distribution of gene expression ( Sahoo et al., 2020 ). SciPy library has been used to calculate kurtosis (k) and skew (s) of the distribution. Supplementary Table Legends Table S1 – List of gene signatures used in ssGSEA, along with their corresponding reference. Table S2 – List of network nodes and edges for the network in Fig 3A , given as an input topology file to RACIPE, along with their corresponding reference. Supplementary Figures Download figure Open in new tab Figure S1: Correlations of RKIP and BACH1 expression with Immune Pathway signatures and Survival in TCGA. Heatmaps representing pairwise Spearman’s correlation coefficient of A) ssGSEA scores of metabolic pathway gene signatures: Glycolysis,) Fatty Acid Oxidation (FAO), Oxidative Phosphorylation (OXPHOS) with each other. B) ssGSEA scores of PD-L1 expression with RKIP and BACH1 C) Heatmap depicting Log2 normalized Hazard Ratios for RKIP and BACH1 across cancers. *: p-value < 0.05. Download figure Open in new tab Figure S2: A) Scatter plots showing correlations of RKIP and BACH1 with SNAI2, VIM and OVOL2 in TCGA breast cancer data. B) Heatmap representing Spearman’s correlation coefficient of ssGSEA scores of RKIP Pathway Metastasis Signature (RPMS) and BACH1 Pathway Metastasis Signature (BPMS) with RKIP, BACH1 and each other, across various cancers. * denotes p-value < 0.05 Download figure Open in new tab Figure S3: A) Scree-plot representing percentage explained variance captured by each principal component from PCA performed on steady state gene expression data obtained from RACIPE. B) PCA Correlation Circle representing correlations of RACIPE steady state gene data along PC1 and PC2 axes. C) i) Kernel Density Estimate plots of z-normalized expression values obtained from RACIPE simulation, for each node in the network ( Figure 3A ). Bimodality coefficients for each distribution are indicated as BC. ii) Pairwise correlation matrix showing Spearman correlation coefficient values of the corresponding pairs of genes. p-value < 10 −4 for all pairwise correlations in the matrix. Download figure Open in new tab Figure S4: Forest plots comparing overall survival (OS) for different combinations of RKIP (low, high) and BACH1 (low, high) in TCGA samples. p-values based on log-rank test, and those with significant differences p < 0.05, p<0.01 and p<0.001 are marked with *, **, and ***, respectively. A) Hazard Ratios (HR) for Overall Survival in MESO (left), SARC (middle) and BRCA (right). B) HR for progression Free Survival in UVM (left), SARC (middle) and BLCA (right). C) HR for Disease-specific Survival in SARC (top) and BLCA (bottom). D) HR for disease-specific Survival in LGG (top) and BRCA (bottom). E) HR for disease-free Survival in THCA (top) and KIRC (bottom). Acknowledgements This work was supported by Ramanujan Fellowship (SB/S2/RJN-049/2018) awarded to MKJ by Science and Engineering Research Board (SERB), Department of Science and Technology, Government of India. Bibliography ↵ Al-Mulla F , Bitar MS , Feng J , Park S , Yeung KC . A New Model for Raf Kinase Inhibitory Protein Induced Chemotherapeutic Resistance . PLoS One 2012 ; 7 : e29532 . https://doi.org/10.1371/JOURNAL.PONE.0029532 . OpenUrl CrossRef PubMed Web of Science ↵ Alsuliman A , Colak D , Al-Harazi O , Fitwi H , Tulbah A , Al-Tweigeri T , et al. Bidirectional crosstalk between PD-L1 expression and epithelial to mesenchymal transition: Significance in claudin-low breast cancer cells . Mol Cancer 2015 ; 14 : 1 – 13 . https://doi.org/10.1186/S12943-015-0421-2/FIGURES/5 . OpenUrl CrossRef PubMed ↵ Beach S , Tang H , Park S , Dhillon AS , Keller ET , Kolch W , et al. Snail is a repressor of RKIP transcription in metastatic prostate cancer cells . Oncogene 2008 27:15 2007 ; 27 : 2243 – 8 . https://doi.org/10.1038/sj.onc.1210860 . OpenUrl PubMed ↵ Brabletz T , Kalluri R , Nieto MA , Weinberg RA . EMT in cancer . Nat Rev Cancer 2018 ; 18 : 128 – 34 . https://doi.org/10.1038/nrc.2017.118 . OpenUrl CrossRef PubMed ↵ Brown MS , Abdollahi B , Wilkins OM , Lu H , Chakraborty P , Ognjenovic NB , et al. Phenotypic heterogeneity driven by plasticity of the intermediate EMT state governs disease progression and metastasis in breast cancer . Sci Adv 2022 ; 8 : eabj8002 . https://doi.org/10.1101/2021.03.17.434993 . OpenUrl ↵ Burns KA , Korach KS . Estrogen receptors and human disease: An update . Arch Toxicol 2012 ; 86 : 1491 – 504 . https://doi.org/10.1007/s00204-012-0868-5 . OpenUrl CrossRef PubMed ↵ Celià-Terrassa T , Bastian C , Liu DD , Ell B , Aiello NM , Wei Y , et al. Hysteresis control of epithelialmesenchymal transition dynamics conveys a distinct program with enhanced metastatic ability . Nat Commun 2018 ; 9 : 5005 . https://doi.org/10.1038/s41467-018-07538-7 . OpenUrl CrossRef ↵ Celià-Terrassa T , Kang Y. Distinctive properties of metastasisinitiating cells . Genes Dev 2016 ; 30 : 892 – 908 . https://doi.org/10.1101/gad.277681.116.892 . OpenUrl Abstract / FREE Full Text ↵ Cessna H , Baritaki S , Zaravinos A , Bonavida B. The Role of RKIP in the Regulation of EMT in the Tumor Microenvironment . Cancers (Basel) 2022 ; 14 : 4596 . https://doi.org/10.3390/CANCERS14194596/S1 . OpenUrl ↵ 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 . https://doi.org/10.3389/fbioe.2020.00220 . OpenUrl ↵ Chauhan L , Ram U , Hari K , Jolly MK . Topological signatures in regulatory network enable phenotypic heterogeneity in small cell lung cancer . Elife 2021 ; 10 : e64522 . https://doi.org/ https://doi.org/10.7554/eLife.64522 . OpenUrl Chen Limo , Gibbons DL , Goswami S , Cortez MA , Ahn YH , Byers LA , et al. Metastasis is regulated via microRNA-200/ZEB1 axis control of tumour cell PD-L1 expression and intratumoral immunosuppression . Nat Commun 2014 ; 5 . https://doi.org/10.1038/ncomms6241 . ↵ Chen Y , Jin L , Ma Y , Liu Y , Zhu Q , Huang Y , et al. BACH1 promotes lung adenocarcinoma cell metastasis through transcriptional activation of ITGA2 . Cancer Sci 2023 . https://doi.org/10.1111/CAS.15884 . ↵ Cho JH , Park S , Kim S , Kang SM , Woo TG , Yoon MH , et al. RKIP Induction Promotes Tumor Differentiation via SOX2 Degradation in NF2-Deficient Conditions . Mol Cancer Res 2022 ; 20 : 412 – 24 . https://doi.org/10.1158/1541-7786.MCR-21-0373 . OpenUrl ↵ Chung VY , Tan TZ , Tan M , Wong MK , Kuay KT , Yang Z , et al. GRHL2-miR-200-ZEB1 maintains the epithelial status of ovarian cancer through transcriptional regulation and histone modification . Sci Rep 2016 ; 6 : 19943 . https://doi.org/10.1038/srep19943 . OpenUrl CrossRef PubMed ↵ Dangi-Garimella S , Yun J , Eves EM , Newman M , Erkeland SJ , Hammond SM , et al. Raf kinase inhibitory protein suppresses a metastasis signalling cascade involving LIN28 and let-7 . EMBO J 2009 ; 28 : 347 – 58 . https://doi.org/10.1038/EMBOJ.2008.294 . OpenUrl Abstract / FREE Full Text ↵ Davudian S , Shajari N , Kazemi T , Mansoori B , Salehi S , Mohammadi A , et al. BACH1 silencing by siRNA inhibits migration of HT-29 colon cancer cells through reduction of metastasis-related genes . Biomedicine and Pharmacotherapy 2016 ; 84 : 191 – 8 . https://doi.org/10.1016/j.biopha.2016.09.021 . OpenUrl ↵ Dong Y , Lin X , Kapoor A , Gu Y , Xu H , Major P , et al. Insights of rkip-derived suppression of prostate cancer . Cancers (Basel) 2021 ; 13 . https://doi.org/10.3390/cancers13246388 . ↵ Dongre A , Rashidian M , Eaton EN , Reinhardt F , Thiru P , Zagorulya M , et al. Direct and indirect regulators of epithelial–mesenchymal transition– mediated immunosuppression in breast carcinomas . Cancer Discov 2021 ; 11 : 1286 – 305 . https://doi.org/10.1158/2159-8290.CD-20-0603 . OpenUrl Abstract / FREE Full Text ↵ Dongre A , Rashidian M , Reinhardt F , Bagnato A , Keckesova Z , Ploegh HL , et al. Epithelial-to-mesenchymal transition contributes to immunosuppression in breast carcinomas . Cancer Res 2017 ; 77 : 3982 – 9 . https://doi.org/10.1158/0008-5472.CAN-16-3292 . OpenUrl Abstract / FREE Full Text ↵ Han W , Zhang Y , Niu C , Guo J , Li J , Wei X , et al. BTB and CNC homology 1 (Bach1) promotes human ovarian cancer cell metastasis by HMGA2-mediated epithelial-mesenchymal transition . Cancer Lett 2019 ; 445 : 45 – 56 . https://doi.org/10.1016/j.canlet.2019.01.003 . OpenUrl ↵ Hari K , Ullanat V , Balasubramanian A , Gopalan A , Jolly MK . Landscape of epithelial mesenchymal plasticity as an emergent property of coordinated teams in regulatory networks . Elife 2022 ; 11 : e76535 . https://doi.org/10.7554/elife.76535 . OpenUrl CrossRef ↵ Huang B , Jia D , Feng J , Levine H , Onuchic JN , Lu M. RACIPE: A computational tool for modeling gene regulatory circuits using randomization . BMC Syst Biol 2018 ; 12 . https://doi.org/10.1186/s12918-018-0594-6 . ↵ Igarashi K , Nishizawa H , Saiki Y , Matsumoto M. The transcription factor BACH1 at the crossroads of cancer biology: From epithelial–mesenchymal transition to ferroptosis . Journal of Biological Chemistry 2021 ; 297 : 101032 . https://doi.org/10.1016/J.JBC.2021.101032 . OpenUrl ↵ Jia D , Li X , Bocci F , Tripathi S , Deng Y , Jolly MK , et al. Quantifying Cancer Epithelial-Mesenchymal Plasticity and its Association with Stemness and Immune Response . J Clin Med 2019 ; 8 : 725 . https://doi.org/10.3390/jcm8050725 . OpenUrl ↵ Jia D , Park JH , Kaur H , Jung KH , Yang S , Tripathi S , et al. Towards decoding the coupled decisionmaking of metabolism and epithelial-to-mesenchymal transition in cancer . Br J Cancer 2021 ; 124 : 1902 – 11 . https://doi.org/10.1038/s41416-021-01385-y . OpenUrl CrossRef ↵ Jiang P , Li F , Liu Z , Hao S , Gao J , Li S. BTB and CNC homology 1 (Bach1) induces lung cancer stem cell phenotypes by stimulating CD44 expression . Respir Res 2021 ; 22 . https://doi.org/10.1186/S12931-021-01918-2 . ↵ Jolly MK , Celia-Terrassa T. Dynamics of Phenotypic Heterogeneity Associated with EMT and Stemness during Cancer Progression . J Clin Med 2019 ; 8 : 1542 . https://doi.org/10.3390/jcm8101542 . OpenUrl ↵ Jolly MK , Huang B , Lu M , Mani SA , Levine H , Ben-Jacob E. Towards elucidating the connection between epithelial-mesenchymal transitions and stemness . J R Soc Interface 2014 ; 11 . https://doi.org/10.1098/rsif.2014.0962 . ↵ Jolly MK , Tripathi SC , Jia D , Mooney SM , Celiktas M , Hanash SM , et al. Stability of the hybrid epithelial/mesenchymal phentoype . Oncotarget 2016 ; 7 : 27067 – 84 . OpenUrl CrossRef PubMed ↵ Karwacki-Neisius V , Göke J , Osorno R , Halbritter F , Ng JH , Weiße AY , et al. Reduced Oct4 expression directs a robust pluripotent state with distinct signaling activity and increased enhancer occupancy by Oct4 and Nanog . Cell Stem Cell 2013 ; 12 : 531 – 45 . https://doi.org/10.1016/j.stem.2013.04.023 . OpenUrl CrossRef PubMed Web of Science ↵ Kim D , Choi B , Ryoo I , Kwak M-K. High NRF2 level mediates cancer stem cell-like properties of aldehyde dehydrogenase (ALDH)-high ovarian cancer cells: inhibitory role of all-trans retinoic acid in ALDH/NRF2 signaling . Cell Death Dis 2018 ; 9 : 896 . https://doi.org/10.1038/s41419-018-0903-4 . OpenUrl ↵ Kong D , Banerjee S , Ahmad A , Li Y , Wang Z , Sethi S , et al. Epithelial to mesenchymal transition is mechanistically linked with stem cell signatures in prostate cancer cells . PLoS One 2010 ; 5 : e12445 . https://doi.org/10.1371/journal.pone.0012445 . OpenUrl CrossRef PubMed ↵ Lee J , Yesilkanal AE , Wynne JP , Frankenberger C , Liu J , Yan J , et al. Effective breast cancer combination therapy targeting BACH1 and mitochondrial metabolism . Nature 2019 ; 568 : 254 – 8 . https://doi.org/10.1038/S41586-019-1005-X . OpenUrl CrossRef PubMed Lee Jiyoung , Lee Jinho , Farquhar KS , Yun J , Frankenberger CA , Bevilacqua E , et al. Network of mutually repressive metastasis regulators can promote cell heterogeneity and metastatic transitions . Proc Natl Acad Sci U S A 2014 ; 111 : E364 – 373 . https://doi.org/10.1073/pnas.1304840111 . OpenUrl Abstract / FREE Full Text ↵ Lee UJ , Frankenberger C , Yun J , Bevilacqua E , Caldas C , Chin SF , et al. A prognostic gene signature for metastasis-free survival of triple negative breast cancer patients . PLoS One 2013 ; 8 . https://doi.org/10.1371/journal.pone.0082125 . ↵ Lignitto L , LeBoeuf SE , Homer H , Jiang S , Askenazi M , Karakousi TR , et al. Nrf2 Activation Promotes Lung Cancer Metastasis by Inhibiting the Degradation of Bach1 . Cell 2019 ; 178 : 316 – 329 .e18. https://doi.org/10.1016/j.cell.2019.06.003 . OpenUrl CrossRef PubMed ↵ Liu L , Matsumoto M , Matsui-Watanabe M , Ochiai K , Callens BKK , Nguyen LC , et al. BACH1 Expression Is Promoted by Tank Binding Kinase 1 (TBK1) in Pancreatic Cancer Cells to Increase Iron and Reduce the Expression of E-Cadherin . Antioxidants 2022 ; 11 . https://doi.org/10.3390/ANTIOX11081460/S1 . ↵ Lourenco AR , Ban Y , Crowley MJ , Lee SB , Ramchandani D , Du W , et al. Differential contributions of pre-And post-EMT tumor cells in breast cancer metastasis . Cancer Res 2020 ; 80 : 163 – 9 . https://doi.org/10.1158/0008-5472.CAN-19-1427/653773/AM/DIFFERENTIAL-CONTRIBUTIONS-OF-PRE-AND-POST-EMT . OpenUrl Abstract / FREE Full Text ↵ Lu T , Xu R , Li Q , Zhao J ying , Peng B , Zhang H , et al. Systematic profiling of ferroptosis gene signatures predicts prognostic factors in esophageal squamous cell carcinoma . Mol Ther Oncolytics 2021 ; 21 : 134 – 43 . https://doi.org/10.1016/j.omto.2021.02.011 . OpenUrl ↵ Lüönd F , Sugiyama N , Bill R , Bornes L , Hager C , Tang F , et al. Distinct contributions of partial and full EMT to breast cancer malignancy . Dev Cell 2021 ; 56 : 3203 – 3221 .e11. https://doi.org/10.1016/J.DEVCEL.2021.11.006 . OpenUrl ↵ Muralidharan S , Sahoo S , Saha A , Chandran S , Majumdar SS , Mandal S , et al. Quantifying the Patterns of Metabolic Plasticity and Heterogeneity along the Epithelial–Hybrid–Mesenchymal Spectrum in Cancer . Biomolecules 2022a ; 12 . https://doi.org/10.3390/biom12020297 . ↵ Muralidharan S , Sehgal M , Soundharya R , Mandal S , Majumdar SS , Yeshwanth M , et al. PD-L1 Activity Is Associated with Partial EMT and Metabolic Reprogramming in Carcinomas . Current Oncology 2022, Vol 29, Pages 8285-8301 2022b ; 29 : 8285 – 301 . https://doi.org/10.3390/CURRONCOL29110654 . OpenUrl ↵ Nagy Á , Munkácsy G , Győrffy B. Pancancer survival analysis of cancer hallmark genes . Sci Rep 2021 ; 11 . https://doi.org/10.1038/s41598-021-84787-5 . ↵ Niwa H , Miyazaki J , Smith AG . Quantitative expression of Oct-3/4 defines differentiation, dedifferentiation or self-renewal of ES cells . Nat Genet 2000 ; 24 : 372 – 6 . https://doi.org/10.1038/74199 . OpenUrl CrossRef PubMed Web of Science ↵ Noman MZ , Janji B , Abdou A , Hasmim M , Terry S , Tan TZ , et al. The immune checkpoint ligand PD-l1 is upregulated in EMT-activated human breast cancer cells by a mechanism involving ZEB-1 and miR-200 . Oncoimmunology 2017 ; 6 : e1263412 . https://doi.org/10.1080/2162402X.2016.1263412 . OpenUrl CrossRef PubMed ↵ Pasani S , Sahoo S , Jolly MK . Hybrid E/M phenotype(s) and stemness: a mechanistic connection embedded in network topology . J Clin Med 2021 ; 10 : 60 . https://doi.org/10.1101/2020.10.18.341271 . OpenUrl ↵ Penas C , Apraiz A , Muñoa I , Arroyo-Berdugo Y , Rasero J , Ezkurra PA , et al. RKIP Regulates Differentiation-Related Features in Melanocytic Cells . Cancers (Basel) 2020 ; 12 : 1 – 22 . https://doi.org/10.3390/CANCERS12061451 . OpenUrl ↵ Roca H , Hernandez J , Weidner S , McEachin RC , Fuller D , Sud S , et al. Transcription Factors OVOL1 and OVOL2 Induce the Mesenchymal to Epithelial Transition in Human Cancer . PLoS One 2013 ; 8 : e76773 . https://doi.org/10.1371/journal.pone.0076773 . OpenUrl CrossRef PubMed ↵ Sahoo S , Mishra A , Kaur H , Hari K , Muralidharan S , Mandal S , et al. A mechanistic model captures the emergence and implications of non-genetic heterogeneity and reversible drug resistance in ER+ breast cancer cells . NAR Cancer 2021 ; 3 : zcab027 . https://doi.org/10.1101/2021.03.14.435359 . OpenUrl ↵ Sahoo S , Singh D , Chakraborty P , Jolly MK . Emergent Properties of the HNF4α-PPARγ Network May Drive Consequent Phenotypic Plasticity in NAFLD . J Clin Med 2020 ; 9 : 870 . https://doi.org/10.3390/jcm9030870 . OpenUrl ↵ Saxena K , Srikrishnan S , Celia-Terrassa T , Jolly MK . OVOL1/2: Drivers of Epithelial Differentiation in Development, Disease, and Reprogramming . Cells Tissues Organs 2022 ; 211 : 183 – 92 . https://doi.org/10.1159/000511383 . OpenUrl ↵ Shi L , Dong B , Li Z , Lu Y , Ouyang T , Li J , et al. Expression of ER-α 36, a novel variant of estrogen receptor α, and resistance to tamoxifen treatment in breast cancer . Journal of Clinical Oncology 2009 ; 27 : 3423 – 9 . https://doi.org/10.1200/JCO.2008.17.2254 . OpenUrl Abstract / FREE Full Text ↵ Shvartsur A , Givechian KB , Garban H , Bonavida B. Overexpression of RKIP and its cross-talk with several regulatory gene products in multiple myeloma . Journal of Experimental and Clinical Cancer Research 2017 ; 36 : 1 – 14 . https://doi.org/10.1186/S13046-017-0535-Z/TABLES/3 . OpenUrl ↵ Spizzo R , Nicoloso MS , Lupini L , Lu Y , Fogarty J , Rossi S , et al. miR-145 participates with TP53 in a death-promoting regulatory loop and targets estrogen receptor-α in human breast cancer cells . Cell Death & Differentiation 2010 17:2 2009 ; 17 : 246 – 54 . https://doi.org/10.1038/cdd.2009.117 . OpenUrl CrossRef PubMed ↵ Subbalakshmi AR , Sahoo S , Mcmullen I , Saxena AN , Venugopal SK , Somarelli JA , et al. KLF4 Induces Mesenchymal-Epithelial Transition (MET) by Suppressing Multiple EMT-Inducing Transcription Factors . Cancers (Basel) 2021 ; 13 : 5135 . https://doi.org/10.3390/cancers13205135 . OpenUrl ↵ Subramanian A , Tamayo P , Mootha VK , Mukherjee S , Ebert BL , Gillette MA , et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles . 2005 . ↵ Tan TZ , Miow QH , Miki Y , Noda T , Mori S , Huang RY , et al. Epithelial-mesenchymal transition spectrum quantification and its efficacy in deciphering survival and drug responses of cancer patients . EMBO Mol Med 2014 ; 6 : 1279 – 93 . https://doi.org/10.15252/emmm.201404208 . OpenUrl Abstract / FREE Full Text ↵ Vilchez Mercedes SA , Bocci F , Ahmed M , Eder I , Zhu N , Levine H , et al. Nrf2 Modulates the Hybrid Epithelial/Mesenchymal Phenotype and Notch Signaling During Collective Cancer Migration . Front Mol Biosci 2022 ; 9 : 279 . https://doi.org/10.3389/FMOLB.2022.807324/BIBTEX . OpenUrl ↵ Viswanathan VS , Ryan MJ , Dhruv HD , Gill S , Eichhoff OM , Seashore-Ludlow B , et al. Dependency of a therapy-resistant state of cancer cells on a lipid peroxidase pathway . Nature 2017 547:7664 2017 ; 547 : 453 – 7 . https://doi.org/10.1038/nature23007 . OpenUrl CrossRef PubMed ↵ Wan Y , Cohen J , Szenk M , Farquhar KS , Coraci D , Krzysztoń R , et al. Nonmonotone invasion landscape by noise-aware control of metastasis activator levels . Nature Chemical Biology 2023 2023 : 1 – 13 . https://doi.org/10.1038/s41589-023-01344-z . ↵ Wang Q , Gun M , Hong X yu . Induced Tamoxifen Resistance is Mediated by Increased Methylation of E-Cadherin in Estrogen Receptor-Expressing Breast Cancer Cells . Sci Rep 2019 ; 9 : 14140 . https://doi.org/10.1038/s41598-019-50749-1 . OpenUrl ↵ Wei C , Cao Yuan , Yang X , Zheng Z , Guan K , Wang Q , et al. Elevated expression of TANK-binding kinase 1 enhances tamoxifen resistance in breast cancer . Proc Natl Acad Sci U S A 2014 ; 111 . https://doi.org/10.1073/PNAS.1316255111/-/DCSUPPLEMENTAL/PNAS.201316255SI.PDF . ↵ Wei X , Guo J , Li Q , Jia Q , Jing Q , Li Y , et al. Bach1 regulates self-renewal and impedes mesendodermal differentiation of human embryonic stem cells . Sci Adv 2019 ; 5 . https://doi.org/10.1126/SCIADV.AAU7887 . ↵ Wenzel SE , Tyurina YY , Zhao J , St. Croix CM , Dar HH , Mao G , et al. PEBP1 Wardens Ferroptosis by Enabling Lipoxygenase Generation of Lipid Death Signals . Cell 2017 ; 171 : 628 – 641 .e26. https://doi.org/10.1016/J.CELL.2017.09.044 . OpenUrl CrossRef PubMed ↵ Wiel C , Le Gal K , Ibrahim MX , Sayin VI , Bergo MO , Jahangir CA , et al. BACH1 Stabilization by Antioxidants Stimulates Lung Cancer Metastasis Article BACH1 Stabilization by Antioxidants Stimulates Lung Cancer Metastasis . Cell 2019 ; 178 : 330 – 345 .e22. https://doi.org/10.1016/j.cell.2019.06.005 . OpenUrl CrossRef PubMed ↵ Xie X , Tian L , Zhao Y , Liu F , Dai S , Gu X , et al. BACH1-induced ferroptosis drives lymphatic metastasis by repressing the biosynthesis of monounsaturated fatty acids . Cell Death & Disease 2023 14:1 2023 ; 14 : 1 – 16 . https://doi.org/10.1038/s41419-023-05571-z . OpenUrl ↵ Yesilkanal AE , Yang D , Valdespino A , Tiwari P , Sabino AU , Nguyen LC , et al. Limited inhibition of multiple nodes in a driver network blocks metastasis . Elife 2021 ; 10 . https://doi.org/10.7554/ELIFE.59696 . ↵ Yun J , Frankenberger C a , Kuo W-L , Boelens MC , Eves EM , Cheng N , et al. Signalling pathway for RKIP and Let-7 regulates and predicts metastatic breast cancer . EMBO J 2011 ; 30 : 4500 – 14 . https://doi.org/10.1038/emboj.2011.312 . OpenUrl Abstract / FREE Full Text ↵ Zhang J , Zhou C , Jiang H , Liang L , Shi W , Zhang Q , et al. ZEB1 induces ER-α promoter hypermethylation and confers antiestrogen resistance in breast cancer . Cell Death Dis 2017 ; 8 . https://doi.org/10.1038/cddis.2017.154 . ↵ Zhou JX , Huang S. Understanding gene circuits at cell-fate branch points for rational cell reprogramming . Trends in Genetics 2011 ; 27 : 55 – 62 . https://doi.org/10.1016/j.tig.2010.11.002 . OpenUrl CrossRef PubMed Web of Science ↵ Zhu GD , Liu F , OuYang S , Zhou R , Jiang FN , Zhang B , et al. BACH1 promotes the progression of human colorectal cancer through BACH1/CXCR4 pathway . Biochem Biophys Res Commun 2018 ; 499 : 120 – 7 . https://doi.org/10.1016/j.bbrc.2018.02.178 . OpenUrl ↵ Zou Y , Ding L , Coleman M , Wang Z. Estrogen receptor-alpha (ER-α) suppresses expression of its variant ER-α36 . FEBS Lett 2009 ; 583 : 1368 – 74 . https://doi.org/10.1016/j.febslet.2009.03.047 . OpenUrl CrossRef PubMed Back to top Previous Next Posted July 08, 2023. Download PDF Supplementary Material Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following A systems-level analysis of the mutually antagonistic roles of RKIP and BACH1 in dynamics of cancer cell plasticity Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share A systems-level analysis of the mutually antagonistic roles of RKIP and BACH1 in dynamics of cancer cell plasticity Sai Shyam , R Soundharya , Manas Sehgal , Mohit Kumar Jolly bioRxiv 2023.07.08.547338; doi: https://doi.org/10.1101/2023.07.08.547338 Share This Article: Copy Citation Tools A systems-level analysis of the mutually antagonistic roles of RKIP and BACH1 in dynamics of cancer cell plasticity Sai Shyam , R Soundharya , Manas Sehgal , Mohit Kumar Jolly bioRxiv 2023.07.08.547338; doi: https://doi.org/10.1101/2023.07.08.547338 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Systems Biology Subject Areas All Articles Animal Behavior and Cognition (8154) Biochemistry (19000) Bioengineering (15136) Bioinformatics (45149) Biophysics (22922) Cancer Biology (20045) Cell Biology (27236) Clinical Trials (138) Developmental Biology (14115) Ecology (21352) Epidemiology (2067) Evolutionary Biology (25746) Genetics (16306) Genomics (23783) Immunology (18967) Microbiology (43104) Molecular Biology (18287) Neuroscience (94805) Paleontology (710) Pathology (3023) Pharmacology and Toxicology (5170) Physiology (8253) Plant Biology (16225) Scientific Communication and Education (2113) Synthetic Biology (4630) Systems Biology (10361) Zoology (2417) window.__CF$cv$params={r:'a4c617bed849d0e9',t:'MTc5MjMxMTMwMA==',u:'01a14e141a6e77448123a4949fab8fe1',ut:'bnfg30wQ.DkQrxlxFMtvNxpjDeooV9zDLGIhI7oT1PI-1792311302-1.2.1.1-6uKt64C5ngJrrToCjiz4b1SNB2A9nR5JvmE7QVj_3gbJhiWJI2gDJVMVWkjuo1.9IG0_KNI03v1NRoptRQMtQhvCWjAbS8Suar9Eo5GPdnw',i:60};(function(){if(!document.body)return;var s=document.createElement('script');s.src='/cdn-cgi/challenge-platform/scripts/precursor/main.js';document.head.appendChild(s);})();