Molecular networks governing responses to targeted therapies in cancer cells are complex dynamic systems that demonstrate nonintuitive behaviors. We applied a novel computational strategy to infer probabilistic causal relationships between network components based on gene expression. We constructed a model comprised of an ensemble of networks using multidimensional data from cell line models of cell-cycle arrest caused by inhibition of MEK1/2. Through simulation of a reverse-engineered Bayesian network model, we generated predictions of G1–S transition. The model identified known components of the cell-cycle machinery, such as CCND1, CCNE2, and CDC25A, as well as revealed novel regulators of G1–S transition, IER2, TRIB1, TRIM27. Experimental validation of model predictions confirmed 10 of 12 predicted genes to have a role in G1–S progression. Further analysis showed that TRIB1 regulated the cyclin D1 promoter via NFκB and AP-1 sites and sensitized cells to TRAIL-induced apoptosis. In clinical specimens of breast cancer, TRIB1 levels correlated with expression of NFκB and its target genes (IL8, CSF2), and TRIB1 copy number and expression were predictive of clinical outcome. Together, our results establish a critical role of TRIB1 in cell cycle and survival that is mediated via the modulation of NFκB signaling. Cancer Res; 77(7); 1575–85. ©2017 AACR.
The RAF–MEK–ERK signaling network is a key module of cellular signal integration and transcriptional regulation. Extracellular stimuli (i.e., growth factors, cytokines, chemokines, integrins, etc.) lead to the activation of this network bringing on changes in the global gene expression and cellular outcomes, such as cell growth, proliferation, migration, and survival. Much research has been focused on the role of MAPK pathway in cancer progression and multiple MEK1/2 inhibitors have been developed and tested in the clinic with limited success as single agents (1, 2). It has been well documented by others and our previous work that in cell line models of breast cancer, inhibition of MEK results in cell-cycle arrest but not cell death (3, 4). We have previously reported the activation of a feedback-loop leading to EGFR-dependent upregulation of PI3K signaling in response to MEK1/2 inhibition (3). This feedback mechanism contributes to resistance to MEK inhibitor–induced apoptosis in these cells.
The MEK–PI3K feedback mechanism highlights the complexities of signal transduction networks that govern nonintuitive responses of cellular systems to pharmacologic inhibitors. Reliable network models would be important to make predictions about network responses to signal transduction inhibitors. This would allow to efficiently identify key network nodes that could be targeted therapeutically, in particular, as part of combination therapies. Various modeling technologies have emerged that capture network complexities in different ways (5). Among those, Bayesian network inference models are particularly promising as they are not only capable of predicting complex network structures but also suggest directionality of interactions between network nodes, which leads to in silico hypotheses generation (6, 7).
To build a model of transcriptional and cellular responses to MEK inhibition, we assessed time-dependent changes in mRNA expression profiles and cell-cycle distribution following MEK inhibition in breast cancer cells. Using a novel Bayesian network inference computational engine (6), ensembles of networks were calculated that revealed novel MEK-dependent regulators of the cell cycle and suggested so far unknown mechanisms of pathway cross-talk with the NFκB network. These model predictions were experimentally validated in cell culture models and demonstrate a role of one of the MEK-regulated genes, TRIB1, in mediating cross-talk with the NFκB network. Further analyses in a large cohort of human breast cancers provided evidence for the existence of a TRIB1–NFκB network in tumors and revealed TRIB1 as a predictor of breast cancer–free survival.
Materials and Methods
The following reagents were used: U0126 (Promega), EGF (Millipore), mimosine (Sigma), rhTRAIL (Millipore), rhTNFα (Life Technologies), TriplePrep Kit (GE Healthcare). ON-TARGET plus SMARTpools siRNAs, NC (noncoding-negative control oligos) and individual oligos constituting the pools were purchased from Dharmacon. RNAiMax and Lipofectamine LTX transfection reagents were from Invitrogen. The following antibodies were used: αR-TRIB1 (Millipore), αR-CCND1, αM-CCNA2, αR-CDC25A, αR-IER2, αR-pCDK2, αR-pIKKa (Santa Cruz Biotechnology), αM-FLIP (Enzo Life Sciences), αM-BID (BD Biosciences), αM-DR5 (R&D Systems), all other antibodies were from Cell Signaling Technology.
The TRIB1–EGFP construct was a generous gift from Dr. Kiss-Toth (University of Sheffield, Sheffield, United Kingdom). The cyclin D1 promoter–containing construct pD1luc WT and mutant promoter constructs D1-κB1/2m, harboring two point mutations in the D1-κB1 (CGCGACCCCC) and the D1-κB2 (CGCGAGTTTT) binding site (introduced point mutations are underlined), were a gift from Dr. Hinz (Max-Delbrück-Center for Molecular Medicine, Berlin, Germany). AP-1–mutant (AP1m) and EtsA/EtsB–double mutant (EtsA/Bm) CCND1 promoter constructs were generated by site-directed mutagenesis of pD1LucWT construct. NFκB-Luc, pMetLuc-C vector, SEAP vector reporter constructs were from Clontech.
NFκB promoter reporter assay
Cells were cotransfected with 500 ng of NFκB -Luc and 250 ng of SEAP transfection-control vectors for 18 hours and then treated with 10 ng/mL TNFα. Activation of NFκB promoter was assayed using Ready-To-Glow Dual Secreted Reporter Assay system (Clontech) according to the manufacturer's instructions 24 hours after TNFα treatment.
MDA-MB-231, SUM149PT, MDA-MB-436, MDA-MB-468 triple-negative basal breast cancer cell lines were obtained from ATCC and authenticated before experimental work began by single tandem repeat analysis at 15 different gene loci and amelogenin (Genetica). Cell line authentication was performed by J.W. Gray and colleagues. Information about cell culture conditions as well as the source, authentication, clinical, and pathologic features of tumors used to derive the breast cancer cell lines used in this study was described in detail previously (8).
siRNA treatment and synchronization
The cells were transfected with 50 nmol/L specified siRNA pools or noncoding control, according to the manufacturer's instructions using RNAiMax (Invitrogen) transfection reagent. Four hours after transfection, the medium was replaced to the one containing 0.4 mmol/L mimosine for 16 hours. Cells were released from blocking and allowed to progress through the cell cycle for 12 hours, after which the cells were reblocked with mimosine for 12 hours. Cells were collected at 10 hours post-mimosine release for cell-cycle analysis. Cell lysates for RNA and protein extraction were collected at 0 and 10 hours postrelease from mimosine block.
Cell cycle, apoptosis analysis, and immunoblots
Cell cycle and apoptosis analysis were performed by FACS as well as standard immunoblots were generated as described previously (3).
Real-time quantitative RT-PCR
Total RNA was extracted from cells at 24 and 72 hours post-siRNA transfection using RNeasy Micro Kit (Qiagen). It was reverse transcribed to cDNA and quantitative RT-PCR analysis using the TaqMan assay (Applied Biosystems) was performed at Genome Analysis Core Facility of Helen Diller Family Comprehensive Cancer Center (UCSF, San Francisco, CA). PCR primers and TaqMan probes for CCND1, TRIB1, IER2, CDKN2C, NUAK1, C14ORF133, CCNE2, TBK1, EGR1, NPC1, SPRED2, KIAA0649, DR5, and YY1 were purchased from Applied Biosystems. hGUS was used as a normalization control. The details of qPCR are described in Supplementary Methods.
Transcriptional assessment of MEK inhibition
We assessed the temporal changes in gene expression profiles induced by EGF and UO126 in the MDAMB231 cell line by RNA expression array technology as described in Supplementary Methods.
TaqMan low-density array
Quantitative PCRs were performed using custom-made TaqMan low-density array (TLDA) focused on 42 cell-cycle–related genes (Applied Biosystems; Supplementary Table S3). Normalization of cDNA levels in the different wells was done using the GAPDH and hGUS RNA as a reference. The amplification protocol and data analysis were carried out as described previously (9). Initial data analysis was performed using Qbase-plus software from Biogazelle (http://www.biogazelle.com/). Each sample set was analyzed using the standardization procedure as described previously (10).
Microarray RNA expression analysis
Cell line and treatment preparation.
Twenty-four hours prior to treatment, cells were transferred to low serum conditions (0.1% FBS), and then they were subjected to the following course of treatment. First, we pretreated the cells with U0126 (10 μmol/L) or DMSO (vehicle control) for 30 minutes (time −0.5 hour). At time 0 hours, the cells were stimulated with EGF (10 ng/mL). Cells were harvested at 8 timepoints: 0, 0.5, 1, 4, 8, 12, 24, and 48 hours after EGF stimulation (two biological replicates). Total RNA was extracted using the RNeasy Mini kit (Qiagen). Microarray data were generated at the LBNL Molecular Profiling Laboratory using a high-throughput, automated GeneChip HT System (Affymetrix). HT_HG-U133A array plates consisting of 96 arrays were used. An Affymetrix liquid handling system (GCAS) was used to process the one-cycle IVT target preparation, array plate hybridization setup, washing, and staining. Plate scanning was performed using a CCD-based high-throughput scanner (Affymetrix). Expression data have been deposited in the NIH Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/; accession number: GSE61364).
Processing of microarray data, gene selection, and data frame construction
The quality of the microarray data was assessed using the Simpleaffy package from R/Bioconductor. Microarray data were then normalized using the Robust Multi-array (RMA) method. A subset of genes was then selected to construct the gene regulatory network based on analysis of covariance model to evaluate the significance of each gene in explaining cell-cycle progression. The 149 top-ranked genes were identified for inclusion. This list was further supplemented by 20 known cell-cycle genes, whose expression was inhibited by treatment with U0126. To learn relationships between drug treatment, mRNA expression profiles, and cell-cycle distribution, gene expression measurements at 0, 0.5, and 24 hours were matched with cell-cycle distribution data at 0, 24, and 48 hours post-treatment.
Learning an ensemble of probabilistic models from data
An ensemble of Bayesian network models including drug treatment, gene expression, and endpoint readouts was learned from the data. Bayesian networks provide a convenient framework to represent the global joint probability distribution of P(X1…….Xn: Q), where X1, … Xn are variables in the model and Θ are parameters. As described previously (11), the first phase of learning proceeds by considering all possible combinations of drug, transcripts, and phenotypic data to obtain a collection of highly likely network fragments. Network fragments representing quantitative relationships among all possible sets of two and three measured variables were considered. Each fragment was assigned a Bayesian probabilistic score indicative of how likely the given fragment was given the data, penalized for its mathematical complexity. Here, the Laplace approximation (11) was used to integrate over posterior parameter distributions to compute the Bayesian scores for each fragment. In the second phase of learning (6), the network fragments were combined to form an ensemble of models. An ensemble is defined as a statistical sample from the distribution P (Model|Data), the probability distribution over all models given the observed data. A network's score was computed as the sum over the scores of all network component fragments. A parallel global optimization procedure employing simulated annealing (the Metropolis–Hastings algorithm) was used to construct an ensemble of 1,024 networks that correctly sampled the probability distribution given the data (6, 12).
Model interventional simulations and statistical analysis
Stochastic simulations of the models were used to generate predictions on the distribution of a variable under different conditions (11). Simulations were performed on the ensembles to predict the distribution of model variables under different perturbations. Then, the predicted distribution of networks in response to the perturbation was compared with the predicted distribution under baseline expression level of the perturbed transcript. A Student t test was used to compute the significance of a transcript's influence on the phenotype and the transcripts were ranked by their P values.
Identification of drivers of cell-cycle progression
To identify genes that were drivers of cell-cycle progression, a systematic perturbation procedure was employed. Genes in the model were set to baseline as well as 10-fold below baseline, and cell-cycle distribution (%G1) was predicted. All genes were evaluated with this procedure and the significance of each gene was evaluated via a Student t test. The genes were then ranked based on their P values, which correspond to their likeliness to be causal drivers of cell-cycle progression. The P value was derived on the basis of 30 random samples from the posterior distribution under each simulation condition.
Gene expression, copy number, and survival analysis of primary tumors
To evaluate the association between TRIB1-related signaling and clinical characteristics, we exploited a large, clinically annotated breast cancer cohort for which paired Illumina HT-12 expression data and Affymetrix SNP 6.0 copy number data were available (13). Here, we report on a subset of 1,980 cases to evaluate the correlation between TRIB1 and NFκB/TRAIL signaling as well as clinical outcome associations. In particular, we computed the Pearson correlation between a panel of 38 genes highlighted by functional analyses of breast cancer cell lines and Bayesian network analysis, and reported to be associated with NFκB/TRAIL signaling. Details of this analysis are provided in Supplementary Methods.
A Bayesian network inference model identifies genes involved in G1 arrest in response to MEK inhibition
To assess the effects of MEK inhibition on global gene expression profiles and cell-cycle distribution, we treated MDA-MB-231 breast cancer cells with U0126, a small-molecule inhibitor of MEK1/2, (Supplementary Fig. S1). To infer relationships between drug treatment, mRNA expression profiles, and cell-cycle distribution, gene expression measurements at 0, 0.5, and 24 hours were matched with cell-cycle distribution data at 0, 24, and 48 hours post-treatment. Initial statistical analysis of the data identified 149 genes that were significantly correlated with cell-cycle arrest following drug treatment (Supplementary Table S1). An additional 20 known cell-cycle genes, whose expression levels were significantly changed by U0126, were also included in the analysis (Supplementary Table S1). These data were used to reverse-engineer an ensemble of 1,024 Bayesian networks (Fig. 1). The inclusion of the 20 known cell-cycle genes impacted by U0126 served as a positive control for testing the models as we expected the model predictions to recover some of the known biology. The ensemble of networks was reverse-engineered from the data in two steps (6, 7, 14, 15). First, a Bayesian score, which takes into account the uncertainty in the parameters and penalizes for complexity, was computed for all possible two- and three-element sub-networks. Second, Metropolis Monte Carlo global optimization was used to infer a statistical sample, or ensemble, of network structures from these sub-networks that was consistent with the data (6). Variance in the ensemble reflected uncertainty in network structure given the information available in the data (6). In this analysis, we retained both the network topology structures as well as the parameterization of the sub-networks to generate quantitative predictions on the effect of gene interventions that account for uncertainty in network topologies (6). Next, systematic in silico simulations of 10-fold knockdown of each of the 169 transcripts were completed to predict their effect on G1–S transition. This was accomplished by setting the value of a gene expression variable to 10% of its original value, and then propagating the effects of the intervention through the network structure, and generating out-of-sample simulated values for the G1 fraction under the knockdown condition (Supplementary Fig. S2). As all networks in the ensemble of models are likely given the data, the entire distribution of predictions generated from the ensemble of networks was analyzed for statistical differences between the knockdown and no-knockdown simulation conditions using a paired t test that samples from the distribution 30 times.
Using this approach, we identified 12 genes as candidate regulators of G1–S transition, based on a significance level of P < 0.01 (Supplementary Table S1). These genes included known regulators of G1–S transition [CCND1 (16, 17), CCNE2 (18), CDC25A (19), CDKN2C (20), NUAK1 (21)] as well as novel candidates (TRIB1, IER2, C14ORF133, TAF11, EBAG9, COQ9, and TRIM27).
Experimental validation of model predictions
To validate the model predictions, we assessed cell-cycle distribution in cells treated with pools of siRNAs (Fig. 2A). We utilized pools of siRNA as they have been demonstrated to result in greater phenotypic penetrance and reduced off-target effects compared with single siRNAs (22). Cells were treated with pools of siRNAs directed against the 12 candidate genes and 5 genes (NPC1, EGR1, SPRED2, BTG3, and KIAA0649) that were predicted to have no effect on G1–S phase transition as negative controls, as well as noncoding control siRNAs (Supplementary Table S2). Target knockdown was confirmed by mRNA analysis and protein levels (where antibodies were available; Supplementary Fig. S3). FACS analysis of cell-cycle distribution confirmed the established role of CCND1, CDC25a, and CCNE2 in G1–S progression, and validated model predictions for novel candidate regulators such as IER2, TRIB1, NUAK1, C14ORF133, TRIM27, TAF11, and EBAG9 (Fig. 2A). However, we saw no effect on G1–S transition with siRNAs directed against CDKN2C and COQ9, while knockdown of KIAA0649, a negative control, resulted in G1 arrest. Effect of individual siRNA duplexes against TRIB1 on cell cycle was also analyzed (Supplementary Fig. S4). Overall, the model predicted the effect of interventions with 82% accuracy, 91% sensitivity, and 67% specificity. To understand the mechanism of G1 arrest, we analyzed the status of proteins known to promote G1–S transition (CCND1, CCNA2, CCNE2, phospho-CDK2, ppRB), as well as inhibitors (p27, p18; Fig. 2B and C). We found that knockdown of genes leading to G1 arrest resulted in at least a 2-fold increase in protein levels of p27, p18, or both, with the exception of TRIM27 and KIAA0649. Simultaneously, knockdown of these genes resulted in inhibition of at least three promoters of G1–S transition (Fig. 2C). Analysis of mRNA expression of 46 genes associated with G1–S progression showed distinct patterns of gene expression after siRNA treatment that confirmed the protein findings for CCND1, CCNA2, and CCNE2 (Supplementary Fig. S5).
TRIB1 regulates CCND1 expression via κB and AP1 sites
TRIB1 was predicted by the model to regulate G1–S phase transition with the highest significance score (P = 2.49E−10). Its siRNA-mediated knockdown resulted in inhibition of cyclin D1 RNA and protein levels (Fig. 2B; Supplementary Fig. S5). To investigate the mechanism of TRIB1 knockdown effect on CCND1 expression, we evaluated its effect on the activity of CCND1 promoter using constructs containing mutations in AP1, NFκB, and Ets transcription factor-binding sites (Fig. 3A and B). In cells treated with NC siRNA, inactivation of either both NFκB (D1-κB1/2M) sites or AP1 site led to the inhibition of CCND1 promoter activity, while inactivation of the two Ets sites had no effect. TRIB1 knockdown resulted in 2-fold downregulation of wild-type CCND1 (D1-WT) promoter activity with further reduction in activity of D1-κB1/2M and AP1 promoter constructs, indicating that TRIB1 signals upstream of both ERK1/2 and NFκB (Fig. 3B).
TRIB1 mediates NFκB-responsive promoter activity and expression of CXCL1, CSF2, and IL8
To confirm the role of TRIB1 in regulation of NFκB activity, we analyzed the effect of TRIB1 knockdown and overexpression on an NFκB-responsive promoter construct. This analysis showed that TRIB1 knockdown inhibited the activity of an NFκB-responsive promoter, while TRIB1 overexpression led to a 3-fold increase in promoter activity (Fig. 3C; Supplementary Fig. S6). The expression of siRNA-resistant TRIB1 rescued the phenotype induced by siRNA-mediated TRIB1 knockdown (Supplementary Fig. S7). Furthermore, evaluation of the expression of NFκB target genes involved in metastasis formation (CSF2, CXCL1) and tumor inflammation (IL8) showed that inhibition of TRIB1 with siRNA resulted in the downregulation of these genes on RNA and protein levels (Fig. 3D; Supplementary Fig. S8). Activation of NFκB is mediated by the inactivation of its inhibitors IκBα, p100, and p105 (23) by several upstream kinases, including IKKα. We reasoned that TRIB1 may act upstream of NFκB. TRIB1 knockdown inhibited phosphorylation and degradation of IκBα as well as production of p50 from its precursor p105 (Fig. 3E; Supplementary Fig. S9). In addition, we found that phosphorylation of IKKα/β (Ser176/180), of IKKα (Thr 23), as well as expression of p100 and its processing with generation of p52 were significantly inhibited by downregulation of TRIB1, possibly indicating regulation of p100 expression by TRIB1 (Fig. 3E; Supplementary Fig. S9). These findings provide mechanistic evidence for regulation of NFκB by TRIB1. No changes in RalA phosphorylation at serine 536 were observed. While IKKα and AKT have been shown to phosphorylate RelA, multiple additional kinases are known to phosphorylate RelA at the same site and might maintain phosphorylation following knockdown of TRIB1 (24, 25).
TRIB1 regulates the activity of PI3K–AKT pathway
One of the potential mechanisms of the cytostatic rather than cytotoxic effect of MEK inhibition is the activation of a feedback loop leading to upregulation of AKT signaling (3). In contrast, TRIB1 knockdown leads to simultaneous inhibition of phosphorylation of ERK1/2 and AKT1 on both T308 and S473 (Fig. 4A). Assessment of downstream targets of AKT kinase revealed that TRIB1 knockdown led to an inhibition of phosphorylation of BCL2, GSK3α/β, and FOXO1A/3A (Fig. 4A). In addition, we observed inhibition of phosphorylation of IKKα-T23, an AKT-specific site, providing further evidence for a putative mechanism of NFκB regulation by TRIB1 through its role as a mediator of PI3K–AKT signaling (Fig. 3E).
TRIB1 inhibition upregulates DR5 levels and sensitizes breast cancer cells to TRAIL-induced apoptosis
We postulated that inhibition of antiapoptotic and induction of proapoptotic signals results in increased cell death upon TRIB1 knockdown alone or in combination with the death receptor agonist TRAIL. Treatment of cells with TRIB1 siRNA alone resulted in an increase in apoptosis in all cell lines tested (Fig. 4B). We also observed a 2- to 3-fold increase in apoptosis in cells treated with TRIB1 siRNA in combination with TRAIL treatment (Fig. 4C). These observations were in accord with the induction of PARP cleavage (Fig. 4D). TRAIL-induced apoptosis occurs through an extrinsic mechanism of direct activation of executioner caspases by receptor-activated caspase-8 (26). In addition, apoptosis can be amplified via recruitment of the intrinsic mitochondrial pathway. Activated caspase-8 cleaves proapoptotic protein BID, leading to its translocation to the mitochondria and resulting in mitochondrial apoptosis activation (26). TRIB1 knockdown resulted in increased caspase-8 cleavage (Supplementary Fig. S10), as well as decreased levels of FLIP and BID, indicating increased cleavage of these proteins and the involvement of both pathways (Fig. 4D). To explain how TRIB1 inhibition leads to the upregulation of TRAIL-induced apoptosis, we assayed both mRNA and protein levels of TRAIL Receptor 2 (DR5). DR5 expression is negatively regulated by NFκB via induction of the transcriptional inhibitor YY1 (27). Here we show that knockdown of TRIB1 leads to a significant upregulation of both, DR5 mRNA and protein levels, as well as downregulation of YY1 mRNA (Fig. 4E and F; Supplementary Fig. S11). These data suggest that TRIB1 knockdown sensitizes cells to TRAIL-induced apoptosis by inhibition of NFκB signaling and upregulation of DR5, as well as by inhibition of AKT prosurvival signaling.
Association between TRIB1/TRAIL/NFκB signaling components in human breast cancer
To establish the significance of TRIB1/TRAIL/NFκB signaling in the context of human disease and to perform exploratory tests of clinical outcome associations, we exploited a large cohort (n = 1,980) of primary breast tumors, with long-term clinical follow-up and which have been profiled at the genomic and transcriptomic levels (13). Given our finding that TRIB1 mediates cyclin D1 expression via the regulation of NFκB and AP1, with TRIB1 knockdown resulting in the downregulation of NFκB targets in vitro, we first sought to examine these relationships in the human system. In agreement with our earlier findings, TRIB1 and NFKB1(p105) expression levels were significantly correlated (ρ = 0.102, P < 0.0001), and TRIB1 was correlated with the expression of the NFκB target genes IL8 (ρ = 0.05, P < 0.05) and CSF2 (ρ = 0.05, P < 0.005). We also comprehensively examined the relationships between TRIB1, TRAIL, and NFκB signaling using this tumor cohort to compare the correlations between a panel of 38 genes believed to be associated with these pathways based on the literature (Fig. 5A). In addition to those listed above, we found that the following gene sets were significantly (P < 0.05) positively correlated with TRIB1: proliferation-associated genes (CCND1, TOP2A, MCM2, MKI67), the antiapoptotic genes (BCL2L1, BCL2L11, CFLAR), the proapoptotic genes (BIK, HRK), components of the death receptor complex (TRAF2, TRAF6), as well as TNFRSF10A/DR4 and IL6. A small number of genes in each of these categories were negatively correlated with TRIB1, including PCNA (proliferation), BCL2LA1 (antiapoptotic) and BAD (proapoptotic), TRADD and MPRIP (death receptor complex), and CASP3 (caspase activity; Fig. 5A). These findings are in agreement with our observation that TRIB1 knockdown sensitizes cells to TRAIL-induced apoptosis, resulting in increased cleavage of both BID and CFLAR in breast cancer cells and supports the role of TRIB1 in cell-cycle progression and survival.
Association between TRIB1/TRAIL/NFκB–associated gene expression profiles and TRIB1 copy number with clinical outcome
We next sought to explore the correlation between TRIB1/TRAIL/NFκB–associated genes and patient survival in a human breast tumor cohort. A comparison of Kaplan–Meier curves for TRIB1 copy number amplified versus neutral cases suggests an association with both breast cancer–specific (BCSS; log-rank P < 0.01) and overall survival (OS; log-rank P < 0.01; Fig. 5B; Supplementary Table S4). To visualize the association between TRIB1 or other single gene expression models and outcome in an exploratory analysis, continuous expression levels were stratified into tertiles, and the predicted survival curves for individuals in the upper (66.66%) and lower (33.33%) groups were compared using Kaplan–Meier plots and the log-rank test. These comparisons suggest an association between both TRIB1 and IL8 expression and outcome, as assessed by both BCSS (log-rank P < 0.01) and OS (log-rank P < 0.05; Supplementary Fig. S12A–S12D).
To investigate these relationships further, while adjusting for clinically relevant covariates, we examined survival as a function of continuous gene expression data for single and multiple gene models based on the panel of 38 TRIB1/TRAIL/NFκB–associated genes using Cox regression. In particular, multivariable analysis was performed by building Cox proportional hazards models, which included the most relevant clinical variables as covariates, namely, age, number of lymph nodes positive, tumor size, and grade, stratified by ER status, and tumor bank to test for differences in BCSS and OS associated with these genes, as described in Materials and Methods. After adjusting for these variables, only BIK, MPRIP, IL8, and TOP2A were significant (Wald test, P < 0.05) in the full cohort, where BCSS was the outcome (Supplementary Fig. S13; Supplementary Table S5). For OS, BIK, MPRIP, and TOP2A were again significant, as were HRK and CASP3, whereas IL8 was not associated with this outcome (Supplementary Fig. S14; Supplementary Table S6).
MEK-dependent signaling is a highly complex process as highlighted by our previously described finding of an EGFR-dependent feedback loop resulting in activation of the PI3K pathway following MEK inhibition (3). Studies broadly assessing protein expression changes following MEK inhibitor treatment corroborate these findings (28).
Among the best documented cellular consequences of MEK inhibition is induction of cell-cycle arrest, predominantly at the G1–S phase transition (3, 29, 30). This cytostatic effect is mediated by a variety of mechanisms (reviewed in ref. 31) including repression of cyclin D1 (32). The work presented here sheds new light onto the intermediary steps downstream of MEK and upstream of CCND1. We utilized an advanced Bayesian network inference engine to deepen our understanding of the molecular networks driving responses of cancer cells to MEK inhibition.
Bayesian network inference models allow for reconstruction of molecular networks and for making causal predictions about the relationships of individual network components. While the Bayesian approach has been applied to the reconstruction of small signal transduction networks before (33, 34), our study involves automated construction of large ensembles of networks and extensive experimental validation. Our model predictions correctly identify the established role of CCND1, CCNE2, CDC25a, and ARK5/NUAK1 in the regulation of G1–S progression. Moreover, we identified and experimentally validated TRIB1, IER2, C14ORF133, TAF11, EBAG9, and TRIM27 as novel regulators of G1–S transition. The high rate of validated predictions underscores the power of simulatable Bayesian models.
In the model, TRIB1 was the highest ranked gene (P = 2.49E−10) in our list of genes predicted to affect cell-cycle distribution. During Drosophila development, Tribbles has been demonstrated to regulate G2–M transition during anlage formation through a cdc25-dependent mechanism (35, 36). TRIB1 has recently been shown to control differentiation of M2-like macrophages in C/EBPα-dependent manner (37), and it has been implicated in the proliferation and chemotaxis of vascular smooth muscle cells via regulation of MAPK activity (38). Most importantly, in the context of the current study, TRIB1 and TRIB2 have been demonstrated to be sufficient to induce leukemia in mice when overexpressed (39). However, the role of TRIB1 in regulation of the mammalian cell cycle remains unknown. In this study, we present evidence that TRIB1 is involved in the regulation of G1–S progression by modulating CCND1 expression in human breast cancer cells. CCND1 is highly overexpressed and is a pivotal driver of G1–S progression in these cells (40). The CCND1 promoter contains binding sites for multiple transcription factors including four Ets, two NFκB, and one AP1 site (41–43). We found that TRIB1 regulates NFκB and AP1-dependent transactivation of the CCND1 gene. These findings support our observation that TRIB1 knockdown leads to the downregulation of CCND1 protein and RNA levels. Indeed, TRIB1 was found to act upstream of the MAPKK SEK1 in C. elegans (44) and it has been shown to physically interact with MEK1 and p65 subunit of NFκB in mammalian cells and to promote the induction of proinflammatory cytokines in lipocytes (45). Moreover, NFκB has been demonstrated to enhance AP-1–dependent transactivation in lymphocytes (46). Thus, our data illuminate a role for TRIB1 in regulating MAPK-dependent signals, and, in a broader sense, in oncogenesis and inflammation. Interestingly, our data demonstrate that TRIB1 RNA expression might be regulated through the MAPK pathway, suggesting a feed-forward mechanism by which signals through the MAPK pathway increase TRIB1 expression which, in turn, further stimulates signaling through the pathway. This possibility is subject of ongoing investigation. Other TRIB family members, in particular TRIB3, have been shown to interact with MAPK and Notch signaling in breast cancer, highlight these pseudo-kinases as modulators of key oncogenic signaling pathways (47).
Recent studies revealed the upregulation of TRIB1 during acute and chronic inflammation in white adipose tissue of mice (48). In agreement with these findings, we observed inhibition of CXCL1, CSF2, and IL8 expression in response to treatment with TRIB1 siRNA, as well as upregulation of these genes in cells overexpressing TRIB1 construct. It remains to be determined whether TRIB1 serves as a nuclear coactivator of p65 subunit in mammary epithelial cells. However, we found that knockdown of TRIB1 leads to the inhibition of pIKKα, pIκBα, and subsequent pIκBα degradation. These findings point to a cytoplasmic role of TRIB1 in regulation of NFκB-dependent transcription. Phosphorylation of IKKα on Tyrosine-23 has been shown to be strongly dependent on AKT (49). In addition, we observed inhibition of multiple components of the PI3K pathway in response to the knockdown of TRIB1, providing evidence that TRIB1 influences signal transduction of the PI3K–NFκB pathways beyond the nucleus. TRIB1 functions, at least in part, as a scaffold protein (50), suggesting that TRIB1 might mediate protein–protein interactions of key regulators of the PI3K and NFκB pathways, resulting in their activation. This possibility is currently under investigation.
Regulation of apoptosis represents one of the main functions of the NFκB pathway. In solid tumors, the pathway predominantly suppresses apoptosis induction through a multitude of mechanisms (reviewed in ref. 51) and has been shown to be a major contributor to resistance to chemotherapeutic agents (52). Suppression of death receptor signaling by NFκB has been implicated as one possible mechanism (53). Our work demonstrates that knockdown of TRIB1 induces death receptor signaling by downregulation of DR5 transcriptional repressor, YY1, resulting in a significant increase in DR5 levels. Regulation of YY1 by NFκB is well documented (54). Despite promising preclinical results, initial clinical experience with death receptor agonists for solid tumor indications have been disappointing (55). On the basis of our findings, we speculate that high expression levels of TRIB1 contribute to the resistance to this class of drugs, and inhibitors of TRIB1 activity could overcome therapeutic resistance.
Discovery of involvement of TRIB1 in regulating proliferation, apoptosis, and cytokine production suggests that high levels of TRIB1 expression might result in a clinically more aggressive tumor phenotype. In agreement with this hypothesis, we found that increased DNA copy numbers of the TRIB1 gene as well as high level of TRIB1 mRNA expression are associated with poor breast cancer survival in a large cohort of breast cancers (13). TRIB1 is located on chromosome 8q24 in immediate vicinity of the oncogene cMyc, which is often amplified in cancer. Furthermore, we found TRIB1 expression to be associated with elements of the NFκB network (e.g., NFKB1, IL8) that we had identified to be downstream of TRIB1 in our in vitro studies. On the basis of the association of TRIB1 with a key survival endpoint in breast cancer and its biological functions, we are currently pursuing studies to evaluate its potential as a target for therapeutic intervention.
In summary, our study highlights that linear signal transduction cascades, such as the RAS–RAF–MEK–ERK module, are embedded in complex molecular networks. We show that the combination of Bayesian network inference models, in silico simulation, and in vitro validation represents a powerful new approach for discerning nonintuitive relationships and generation of testable hypotheses about signaling networks from a large, multidimensional dataset. By applying this strategy, we identified TRIB1 as an important mediator of key signal transduction pathways and a candidate therapeutic target in breast cancer.
Disclosure of Potential Conflicts of Interest
P. McDonagh has ownership interest (including patents) in a patent. I. Khalil has ownership interest (including patents) in GNS Healthcare. W.M. Korn is a consultant/advisory board member for Merrimack Pharmaceuticals. No potential conflicts of interest were disclosed by the other authors.
Conception and design: R. Gendelman, H. Xing, J.W. Gray, I. Khalil, W.M. Korn
Development of methodology: R. Gendelman, H. Xing, H. Feiler, P. McDonagh, I. Khalil
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): R. Gendelman, O.K. Mirzoeva, P. Sarde
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): R. Gendelman, H. Xing, O.K. Mirzoeva, P. Sarde, C. Curtis, P. McDonagh, W.M. Korn
Writing, review, and/or revision of the manuscript: R. Gendelman, O.K. Mirzoeva, C. Curtis, P. McDonagh, W.M. Korn
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): R. Gendelman, O.K. Mirzoeva, P. McDonagh, J.W. Gray, W.M. Korn
Study supervision: R. Gendelman, I. Khalil, W.M. Korn
We thank E. Kiss-Toth for generating TRIB1-EGFP vector, M. Hinz for providing D1-WT and D1-kB1/2m constructs, B. Church, K. Runge, R. Miller, and D. Decaprio for the development of REFS(TM) algorithms, software, and insights on methodology utilized in this work. We also thank P. Jensen and Boris Hayete for critically reading the manuscript.
This work was supported by NIH, National Cancer Institute grant U54 CA 112970 (J.W. Gray and W.M. Korn) and P30CA082103.
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.