A system-wide approach to monitor responses to synergistic BRAF and EGFR inhibition in colorectal cancer cells

Intrinsic and/or acquired resistance represents one of the great challenges in targeted cancer therapy. A deeper understanding of the molecular biology of cancer has resulted in more efficient strategies, where one or multiple drugs are adopted in novel therapies to tackle resistance. This beneficial effect of using combination treatments has also been observed in colorectal cancer patients harboring the BRAF(V600E) mutation, whereby dual inhibition of BRAF(V600E) and EGFR increases antitumor activity. Notwithstanding this success, it is not clear whether this combination treatment is the only or most effective treatment to block intrinsic resistance to BRAF inhibitors. Here, we investigate molecular responses upon single and multi-target treatments, over time, using BRAF(V600E) mutant colorectal cancer cells as a model system. Through integration of transcriptomic, proteomic and phosphoproteomics data we obtain a comprehensive overview, revealing both known and novel responses. We primarily observe widespread upregulation of receptors tyrosine kinases and metabolic pathways upon BRAF inhibition. These findings point to mechanisms by which the drug-treated cells switch energy sources and enter a quiescent-like state as a defensive response, while additionally reactivating the MAPK pathway.


Introduction
Despite the progressive development of novel drugs for personalized medicine, intrinsic and/or acquired resistance remains a major limitation of targeted anticancer therapies (Groenendijk & Bernards, 2014;Ahronian & Corcoran, 2017). Most of these drugs target components of the mitogen-activated protein kinase (MAPK) signaling pathway, which contains oncogenes such as KRAS, BRAF and the epidermal growth factor (EGFR) (Dhillon et al, 2007;Miyamoto et al, 2017). The use of monotherapy to inhibit these oncogenes has often been found to be ineffective due to reactivation of signaling pathways. For instance, upregulation of upstream components such as receptors tyrosine kinases (RTKs) in KRAS mutant lung and colorectal cancer (CRC) or of downstream components such as KRAS wild-type in CRC have been revealed to be responsible for intrinsic drug resistance (Sun et al, 2014;Karapetis et al, 2008).
Multiple independent studies on CRC found a crucial role of EGFR as a key driver of resistance to BRAFi monotherapy (Prahallad et al, 2012;Corcoran et al, 2012;Klinger et al, 2013). In congruence with the role of EGFR in conferring resistance to BRAFi, the suppression of tyrosine phosphatase non-receptor type 11 (PTPN11)which is required to transduce signals from EGFR and other RTKs to the downstream MAPK pathwayalso sensitizes BRAF(V600E) CRC cells to BRAF inhibition (Prahallad et al, 2015). Consequently, the identification of EGFR as a mediator of intrinsic resistance to BRAFi in CRC has led to initiation of several clinical trials which combine inhibition of both EGFR and BRAF (BRAFi+EGFRi), or of other MAPK pathway members (Sundar et al, 2017).
Although the BRAFi+EGFRi combination treatment is more effective than BRAFi monotherapy in CRC , it remains unclear whether EGFR is the only or most potent synthetic lethal co-target of BRAF(V600E) in CRC. Addressing this issue requires an understanding of the cellular response to drug treatment across different molecular levels. Such multilayer approaches could elucidate different branches of the signaling network, and track how perturbations propagate to gene and protein expression in driving resistance. Several studies have already highlighted the widespread responses to drug treatment in cancer using multi-omics approaches and adequate data integration (Mertins et al, 2016;Zhang et al, 2016Zhang et al, , 2014. Advances in next-generation sequencing and proteomics approaches (Altelaar et al, 2012;Mathivanan et al, 2008) in combination with enhanced data integration solutions have paved the way for such important investigations (Nesvizhskii, 2014;. Notably, the integrated use of transcriptomics and (phospho)proteomics has recently demonstrated its power in describing physiopathological processes through phenotype and proteotype analysis (Faulkner et al, 2015;Low et al, 2013;Cutillas, 2015;Mertins et al, 2016).
In this study we analyze and integrate proteomics, phosphoproteomics, and transcriptomics data to follow the molecular responses over time upon perturbation with BRAFi, EGFRi, or their combination in CRC cell lines (Figure 1). We aim to study whether there are other posttranslational or transcriptional mechanismsin addition to EGFR and the MAPK pathwaythat are activated upon treatment, in order to identify novel targets that may overcome innate and acquired resistance to BRAF inhibition.

Results
A multi-omics overview of BRAF mutated colorectal cancer cell response to targeted drug treatment The WiDr CRC cell line harboring the BRAF(V600E) mutation was selected as a model system for our analyses (Noguchi et al, 1979;Chen et al, 1987). To study the differences in signaling we treated the WiDr cells with either BRAFi or EGFRi monotherapies or with the combination treatment BRAFi+EGFRi. Additionally, we employed a WiDr PTPN11 knockout (KO) cell line treated with BRAFi (BRAFi in PTPN11 KO) to investigate if there are functional differences between PTPN11 KO and EGFRi. Cell growth was synchronized by serum starvation for 24 hours (h), followed by 30 minutes (min) incubation with or without drugs before serum stimulation (Supplementary Methods). Unstimulated control and PTPN11 KO control cells were immediately harvestedindicated as T=0 h throughout this studywhereas stimulated samples were collected in a time course at 2, 6, 24, and 48 h after treatment ( Figure 1B). We performed transcriptomic (RNA-seq) and (phospho)proteomic profiling at each of the time points as indicated in Figure 1C. This study was designed to capture the initial responses to the different drug treatments and the elicited RTK signaling, monitoring their effects on gene and protein expression, and the onset of feedback mechanisms. For this approach, a customized protein sequence database was generated using RNA-seq data to account for WiDr-specific non-synonymous variants ( Figure 1C and Supplementary Methods).
The phosphorylation profile of ERK (MAPK1) (Figure 2A, top panel) was used as a positive control to verify the drug-induced regulation and overall quality of the label-free quantitative (phospho)proteomics approach. In concordance with previous studies (Prahallad et al, 2015), pERK becomes downregulated upon BRAFi, and this effect is enhanced by the addition of EGFRi, reaching, as expected, the same levels as of the BRAFi in PTPN11 KO cell line.
Complementary Western blots (Figure 2A, bottom panel) show excellent correlation with the label-free phosphoproteomics quantified data. Further quality analysis demonstrates high correlation between respective biological replicates for each omics dataset, with median correlations of 0.99, 0.93 and 0.83 for the transcriptomics, proteomics and phosphoproteomics data respectively ( Figure 2B). Quantified proteomics data points show, as expected, slightly higher variability in comparison to RNA-seq data (Haider & Pal, 2013), while phosphoproteomics data points exhibit even higher variability. The final (phospho)proteomics dataset consisted of 5692 quantified protein groups and 7141 quantified Class I phosphosites (localization probability > 0.75), both of which were measured in at least one timepoint of any of the six applied conditions ( Figure S1). The transcriptome dataset contained a final list of 21,446 genes. We developed a Graphical User Interface (GUI) to facilitate rapid data comparison. The GUI enables selection of a specific gene to immediately visualize a comparison of its expression profile at transcriptomic and (phospho)proteomic levels ( Figure 2C).
To extract a more general overview of the effect of the different experiments on the cells, we performed Principal Component Analyses on each omics data type ( Figure 3A). The trend along the first two principal components is similar for all data types. Principal Component 1 (PC1) represents the variation in the measurements over time in the BRAFi treated samples (BRAFi,BRAFi+EGFRi,and BRAFi in PTPN11 KO). This variation is slightly greater in BRAFi+EGFRi and BRAFi in PTPN11 KO samples compared to the samples treated with BRAFi only. Notably, the onset of the variation in the direction of PC1 occurs earlier in the transcriptomic data (after 6 h) than in the proteomics data (after 24 h) reflecting the delay from translation to transcription. Principal Component 2 (PC2) reflects the variation in the measurements of the analysed cells in non-BRAFi treated samples (control, EGFRi and PTPN11 KO control). In the proteomics data, PC2 also clearly captures the variation induced by the PTPN11 KO.

PTPN11 knockout induces post-transcriptional downregulation
To further investigate the segregation of the isogenic cell line PTPN11 KO from PTPN11 wildtype (WT) as observed predominantly at the protein level, we performed differential protein and mRNA expression analyses using a linear model (limma), where the PTPN11 status (KO and WT), the drug treatments (BRAFi, EGFRi and BRAFi+EGFRi) and the time-points (0, 2, 6, 24 and 48 h) were set as coefficients. Our aim was to elucidate how the PTPN11 status affects proteins and mRNAs differently by contrasting differential expression of proteins and mRNAs (Table S1).
At the protein level, we observe the strongest difference in creatine kinases brain-type (CKB), with a log 2 -fold-change of -6.3 in PTPN11 KO cell line compared to PTPN11 WT. Remarkably, CKB shows minimal difference on average in mRNA expression (log 2 -fold-change = -0.65) ( Figure S2A). Interestingly, CKB has been recently considered responsible for promoting survival in WiDr CRC cells, by regulating the reservoir of high-energy phosphates in tissues to sustain intracellular energetic requirements (Loo et al, 2015). Further analysis identified 64 additional proteins that have significant differential protein expression (FDR = 0.1) between PTPN11 KO and PTPN11 WT, a large negative log 2 -fold-change of < -1 at the protein level, and a small log 2 -fold-change at the mRNA level (< protein log 2 fold-change/2) ( Figure 3B).
These 65 proteins also demonstrate a consistent time-course profile in which expression at both the protein level and the mRNA level is initially upregulated in all six experimental conditions ( Figure 3C). At T > 24 h, downregulation of both the mRNA and protein level occurs in BRAFi treated conditions, whereas expression levels in other conditions are upregulated or remain unchanged. This observation suggests that these 65 proteins can be regulated at the mRNA expression level, despite downregulation in an mRNA-expression independent manner upon PTPN11 KO. The consistency in expression profiles across the 65 proteins also suggests that these genes are to some extent co-regulated.
We next investigated whether these 65 proteins are functionally related by performing an enrichment analysis using the MSigDB Hallmarks gene-sets (Liberzon et al, 2015). Our analysis reveals a strong enrichment of the interferon alpha (IFN-α) and gamma (IFN-γ) response genesets (enrichment > 15-fold, p < 10 -9 , Table S1) which are known to suppress cell viability through the JAK/STAT pathway (Ivashkiv & Donlin, 2013). Interestingly, PTPN11 negatively regulates the INF-induced JAK/STAT pathway by dephosphorylating STAT1 on both residue Y701 and S727 (Du et al, 2005). In line with these observation, STAT1 is also significantly downregulated at the protein level in the PTPN11 KO cells (p < 10 -12 ) in our dataset and shows a time-course profile similar to the aforementioned 65 genes ( Figure S2B).

System-wide propagation of drug perturbation
A key aspect of our study was to provide a system-wide understanding of the propagation of cellular responses from signaling (phosphoproteomics) to gene transcription (transcriptomics) and then translation and protein expression (proteomics) in response to BRAF and/or EGFR inhibition. We therefore performed correlation-based hierarchical clustering on the 1500 phosphosites, mRNAs and proteins exhibiting the highest variance within each dataset, which resulted in eight clusters for each omics dataset ( Figure S3). We observe that in the transcript-, phosphosite-and protein clusters, all BRAFi treated samples (BRAFi, BRAFi+EGFRi and BRAFi in PTPN11 KO) exhibit similar clustering profiles distinct from non-BRAFi treated samples.
Next, we investigated the biological function of all the clusters by performing an enrichment analysis based on transcription factor-target gene (transcriptomic clusters) and kinase-substrate (phosphoproteomics cluster) relationships. We also conducted enrichment analyses using the hallmarks gene-sets from MSigBDg (Liberzon et al, 2015) (Table S2) MAPK, PI3K-AKT and NF-κB signaling are involved in the early response The early response phosphoproteomics cluster is enriched for substrates of kinases belonging to the MAPK and PI3K-AKT pathways ( Figure 4A, right panel). An immediate decrease in phosphorylation of MAPK1 (ERK2) and MAPK3 (ERK1) substrates occurs upon BRAF inhibition, which is consistent with dephosphorylation of Y187 and T185 phosphosites on MAPK1 and Y204 on MAPK3 (Figure 2A and S4A). The early response cluster also includes p70S6K (RPS6KB1) substrates. Similarly, this is consistent with the dephosphorylation of the kinase itself (at the S427 phosphosite) after 2-6 h upon BRAFi and BRAFi+EGFRi ( Figure   S4B). PKBα (AKT1) and PKBβ (AKT2) substrates are also enriched in this cluster exhibiting an unexpected deactivation of AKT upon BRAFi monotherapy with respect to what was previously shown (Prahallad et al, 2012). The extent of deactivation could not directly be resolved since we did not detect the relevant phosphosites of AKT1 and AKT2 in our dataset. Therefore, we Surprisingly, NF-κB, which is not typically considered to be downstream of the MAPK pathway (Moynagh, 2005), is amongst the transcription factors most enhanced in the early response cluster. The phosphorylation profile of S907 on NFKB1 is part of the early response cluster ( Figure 4C). NFKB1 S907 phosphorylation induces activity (Cartwright et al, 2016), suggesting that gene expression of the early response cluster is partially driven by inactivation of NF-κB.
Our proteomics data show no substantial variation in the protein expression of NFKB1 under any experimental condition (Figure S4C), indicating that activity of NFKB1 is likely inhibited by dephosphorylation. The residue S907 of NFKB1 is a possible substrate of GSK3β kinase (Demarchi et al, 2003), which is downstream of MAPK in the pathway, and GSK3β substrates are enriched in the early response phosphoproteomics cluster (albeit not significantly after multiple testing correction). This observation suggests that deactivation of NF-κB is a consequence of the observed MAPK pathway deactivation.
Drug treatment affects the cell cycle and cell proliferation at the late cellular response The late response cluster reflects the effect of drug treatment on cell cycle and proliferation, both at the transcriptomics and phosphoproteomics level. The phosphoproteomics late response cluster is highly enriched for substrates of the cell cycle regulators CDK1 and CDK2 ( Figure 4A, right panel). This finding is corroborated by decreased CDK1 ( Figure S4D) and CDK2 ( Figure S4E) protein expression. Similarly, the transcriptomic late response cluster is enriched for target genes of the cell cycle regulators E2F1, E2F4, and E2F3 ( Figure 4B, right panel). The connection between late response phospho-signaling and gene expression is mediated by downregulation of the phospho-residue T821 on RB1 ( Figure 4D), a substrate of CDK2, while total RB1 protein expression remains relatively constant ( Figure S4F). This downregulation induces binding of RB1 to E2F1 thereby inhibiting E2F1 activity (Lentine et al, 2012), and promoting cell cycle arrest (Henley & Dick, 2012) under BRAFi conditions at the late timepoints. Further enrichment analysis of the MSigDB hallmarks gene-sets revealed enrichment for proliferation-related gene sets including E2F targets, genes involed in the G2M checkpoint, and mitotic spindle genes (Table S2).
We also observe a late response cluster in the proteomics data, which is strongly enriched for targets of MYC and E2F (proteomics cluster 4). MYC and E2F targets are also enriched in the transcriptome early response cluster, again indicative of the delay between transcription and translation.
Combined BRAFi and EGFRi treatment leads to a distinct metabolic response In addition to the early and late response clusters, we observe a set of transcripts, proteins and phosphosite clusters that are upregulated in BRAFi treated samples compared to control and EGFRi-only treated samples. These clusters include phosphoproteomics cluster 5, transcriptomics clusters 5 and 6, and proteomics cluster 2 ( Figure S3) and exhibit the strongest regulation at 48 h. mRNA and protein clusters in this set were enriched for metabolic processes including peroxisome, fatty acid and bile acid metabolism (Table S2). Notably, mRNA cluster 6 is enriched for targets of the transcription factors C/EBPβ, Sp1, HNF-1α and HNF-4α, which are involved in fatty acid metabolism and in the gluconeogenesis pathway (Desvergne et al, 2006).
To pinpoint in more detail which biological processes are upregulated in the BRAFi+EGFRi treated samples at 48 h, we performed differential expression and enrichment analysis of transcripts and proteins that were significantly upregulated at this time point with respect to control (log 2 -fold-change > 1, FDR < 0.05). KEGG and Reactome pathway analyses revealed significant upregulation in metabolic processes and pathways including the pentose phosphate pathway (PPP), the TCA cycle and the lipid metabolism (Table S3). Taken together, these observations suggest that the cells reprogram to a distinctive metabolic pattern in response to the treatment.
Similarly, cluster 5 of phosphoproteomics data was enriched in substrates of the calcium/calmodulin dependent protein kinase II isoforms (CAMKIIα/β/γ/ δ) and of pyruvate dehydrogenase kinase isozyme 1 (PDHK1). Besides being responsible for Ca 2+ homeostasis under pathophysiological conditions (Anderson, 2011), CAMKII is also linked to resistance to apoptosis due to its metabolic activation caused by increased levels of acetil-CoA and intermediates in the PPP (McCoy et al, 2013;Huang et al, 2014). PDHK1 rather regulates glucose and mitochondrial metabolism by phosphorylating pyruvate dehydrogenase E1 component subunit alpha (PDHA1) on serine residues (Harris et al, 2002;Dupuy et al, 2015). In our data, PDHA1 expression is constant at the protein level in the BRAFi+EGFRi samples at 48 h ( Figure S5C), while phosphorylation is upregulated at S232 on PDHA1 ( Figure S5A) but downregulated at S293 ( Figure S5B). As PDHK1 is not detected at the protein level but only at mRNA level ( Figure  Feedback responses aim to reactivate the MAPK pathway Cells are able to preserve a stable homeostatic balance through activation of feedback mechanisms in a counteractive manner as a response to perturbation. Here, we sought to elucidate the molecular response to drug treatment by investigating the presence of potential feedback response mechanisms. In general, a homeostatic response to pathway stimulation can be either downregulation of a signal transducer or upregulation of a signal inhibitor. To study this systematically, Legewie et al. collected gene expression profiles of responses to pathway stimulation (Legewie et al, 2008). By plotting the expression changes of known signal transducers and inhibitors of major signaling pathways (MAPK, PI3K, cAMP, TGFβ, JAK/STAT) against their mRNA half-lives, they observed a striking pattern. All responding genes where a) signal inhibitors and b) short-lived. They called these short-lived signal inhibitors Rapid Feedback Inhibitors (RFIs), which, according to them, constitute a fast and efficient homeostatic response mechanism. We here investigated if a similar design principle applies to pathway inhibition, plotting the log 2 -fold change mRNA expression after 2 hcompared to control T=0 h in each conditionagainst the mRNA half-lives obtained by Legewie et al. (Figure 5A). We  (Figure S7). ERBB2 and ERBB3 are known to confer acquired resistance mechanisms as evidenced by reduced response to EGFR and BRAF(V600E) inhibitors (Montero-Conde et al, 2013;Sun et al, 2014;Jain et al, 2010;Temraz et al, 2016). Interestingly, in addition to the ERBB2/ERBB3 upregulation, the regulator ERBB receptor feedback inhibitor 1 (ERRFI1), which interferes with ERBB family member homo-and hetero dimer formation (Liu et al, 2012), is downregulated in all BRAFi treated samples ( Figure 5E). Altogether, our data suggests the existence of an additional mechanism through which WiDr CRC cells may activate ERBB signaling to compensate for MAPK pathway inhibition.
Combined BRAFi and EGFRi treatment induces a stable but reversible growth arrest To further investigate the response of BRAF(V600E) mutant cells to combined BRAFi+EGFRi treatment, we exposed WiDr cells to the combination of both drugs for a prolonged period of 78 days, after which few surviving cells were still present in small colonies that did not further expand, indicative of growth arrest ( Figure 6A, panel 1). To elucidate whether the arrested cells retained the ability to resume their cell cycle upon drug withdrawal, we switched the cells to conventional medium ("drug off"). As shown in Figure 6A

ERBB inhibitors provide limited benefit in BRAF(V600E) CRC treatment
Next we sought to determine if the observed upregulation of ERBB2 and ERBB3 upon BRAFi+EGFRi treatment could be further exploited. We first studied whether inhibition of ERBB2 and ERBB3 in combination with BRAFi and EGFRi may lead to complete cell death by using gefitinib, lapatinib and sapitinib as known tyrosine-kinase inhibitors of EGFR, EGFR/ERBB2 and EGFR/ERBB2/ERBB3, respectively. Next, we explored if there was an optimal drug concentration, which would be synergistic with the maximum tolerated dose of PLX4032 (BRAFi) (Supplementary Methods). As expected (Prahallad et al, 2012), WiDr CRC cells are resistant to monotherapy of either gefitinib or lapatinib or sapitinib, with decreasing viability only at very high, i.e. cytotoxic, concentrations ( Figure 7A). All three viability curves depict a 60 % decrease in cell viability upon addition of 3 µM PLX4032, with limited benefit from the combination with EGFRi. We did not observe a significant difference in growth inhibition across the different double treatments, suggesting that additional inhibition of ERBB2 and ERBB3 does not provide further synergy with BRAFi.
Inhibitors of metabolic enzymes provide limited benefit in BRAF(V600E)

CRC treatment
We further evaluated combination treatments to target MAPK pathway together with the TCA cycle or with the fatty acid β-oxidation. For this porpuse we used two readily available metabolic drugs: etomoxir, a CPT1 inhibitor, and dichloroacetate (DCA), a pan-inhibitor of PDKs, to determine if inhibition of CPT1 or PDKs enhances sensitivity to therapy with PLX4032 and gefitinib (EGFRi). CPT1 is directly upstream of CPT2 in the peroxisomal fatty acid oxidation processes and is responsible for transporting long chain fatty acids from the cytosol to the mitochondrial matrix (Lodhi & Semenkovich, 2014;Pucci et al, 2016). CPT2 is significantly upregulated in our experiments, nevertheless, to the best of our knowledge, a direct inhibitor of CPT2 has not been previously reported. PDKs are responsible for deactivation of PDHA1 through phosphorylation of serine residues in PDHA1 (Harris et al, 2002;Dupuy et al, 2015) ( Figure S5). We hypothesized that inhibition of PDKs can increase the mitochondrial oxidative state and consequently the amount of reactive oxidative species (ROS) in the cytoplasm causing apoptosis due to high toxicity. We therefore evaluated the combination treatment of each metabolic inhibitor (etomoxir and DCA) with BRAFi and EGFRi on WiDr cell viability. After assessing IC50 of both metabolic drugs (Figure S8A), we selected the clinical doses of 10 µM etomoxir (Ito et al, 2012;Holubarsch et al, 2007;Samudio et al, 2010) or 1 mM DCA (Dunbar et al, 2014;Fox et al, 1996) and we performed two dose-response curves in the presence of 3 µM gefitinib and increasing concentrations of PLX4032. The first curve was obtained by adding the metabolic drug simultaneously to the BRAFi and EGFRi (T=0 h) (Figure S8B), and the second by adding it 96 h after BRAFi and EGFRi (T=96 h) (Figure 7B), when metabolism is supposed to be significantly upregulated according to our data. In both cases, we do not observe any significant differences in the viability curves of the triple treatments in comparison to the double treatments. These findings suggest that inhibition of CPT1 or PDKs does not increase sensitivity to therapy with PLX4032 and gefitinib.

Discussion
In this study we performed an integrative quantitative multi-omics analysis to obtain a systemwide molecular characterization of signaling perturbation over time in WiDr CRC and WiDr PTPN11 KO cell lines, after drug inhibition targeting either BRAF(V600E) and/or of EGFR. Our data reveal that all samples treated with BRAFi show similar response in both PTPN11 WT and KO, with a less pronounced effect following BRAFi-monotherapy. This indicates that main signalling responses depend on the inhibition of BRAF(V600E), whereby the additional inhibition of EGFR by drug further amplifies the effect. Additionally, EGFRi-only treated cells exhibit similar responses to PTPN11 KO samples, confirming that suppression of this secondary signaling pathway confers sensitivity to BRAFi in CRC (Prahallad et al, 2015).
By comparing proteomics and transcriptomics data, we identified a set of genes that are exclusively downregulated at the protein level, upon PTPN11 KO. These proteins are negative regulators of the interferon pathway (Porritt & Hertzog, 2015), involved in controlling immune response. Downregulation of negative regulators may support the immune response elicited by PTPN11 in vivo (Mainardi et al., unpublished data). This finding might be relevant for the development of therapeutic approaches that aim to inhibit PTPN11 activity.
The integrated omics analysis enabled us to trace the system-wide drug response upon treatment, from signalinge.g. inactivation of kinases downstream of the MAPK pathwaythrough transcriptione.g. inhibition of genes downstream in the MAPK pathway. Shutting down MAPK signaling results in downregulation of CDK signaling, inducing cell cycle arrest at a later stage. Besides the inactivation of the MAPK pathway, all three datasets show an increase of oxidative metabolic processes, with significant upregulation of enzymes involved in lipid metabolism and the TCA cycle. We further observed that treatment with BRAFi induces upregulation of RTKs, including ERBB2 and ERBB3, which was found to be more pronounced when co-treated with EGFRi or in PTPN11 KO cells.
A subset of cells survived the BRAFi+EGFRi treatment and could start to proliferate again after drug removal. This indicates that these cells survive by utilizing different metabolic regimes, pointing at potential future avenues on how to target these cells. However, in our work combining inhibition of the MAPK pathway and specific metabolic processes did not result in any significant difference in cell viability. The precise protective role of metabolic adaptation in the ability of cells to tolerate drug treatment remains elusive, and further studies are required.
Apart from the upregulation of metabolic processes, all adaptive responses we observed appear to be homeostatic responses that try to re-activate the MAPK pathway, but their attempting seem to fail in the range time of this study. Despite apparent activation of ERBB family members, we do not find an additional benefit of inhibiting ERBB2 and/or ERBB3 in combination with EGFRi and BRAFi, suggesting that these homeostatic responses are not necessarily functional under the tested conditions. Further studies are required to establish the more general implications of these findings, investigating different cell lines and several medium conditions that more closely mimic physiological environments. Importantly, we do not find any evidence of parallel signaling pathways being activated in response to drug treatment.
Together, this suggests that reactivation of the MAPK pathway is a requirement for BRAF(V600E) mutation CRC cells to become resistant to BRAF inhibition. This view is supported by observations in the clinic that resistance to MAPK pathway inhibitors is typically mediated by mutations or amplification in the MAPK pathway in CRC patients (Ahronian et al, 2015).
Our findings on metabolic rewiring do not show any direct impact when targeted, on the time scale measured in this study, but might still be relevant to therapy as two recent studies in melanoma demonstrated the dependence of resistant cells on mitochondrial respiration (Corazao-Rozas et al, 2013;Cesi et al, 2017). These studies highlight the importance of studying the complete omics landscape to identify relevant drug targets. The integrative multi omics approach employed here provides a time based and in depth view of the signaling mechanisms involved in drug response. Our findings highlight the importance of measuring these different levels simultaneously as exemplified by the RTKs regulation and PTPN11 specific signals. We expect to enable and accelerate future research into these mechanisms by making the data resource available.

Materials and Methods
Experimental design

(Phospho)proteomics
Cell lysates were first digested, then desalted and finally analyzed via mass spectrometry.

Differential expression analysis
All differential expression analyses were performed using limma (Ritchie et al, 2015). mRNA read counts were first transformed using voom (Gentleman et al, 2004). For the comparison of

Clustering
Hierarchical clustering was performed on the 1500 phosphosites, proteins, or mRNAs with the highest variance over all conditions. The pairwise distance was calculated based on the Pearson correlation. The obtained distances were used for hierarchical clustering and the resulting tree was split into 8 groups for each dataset. Enrichment analysis were done using Fisher's exact test. KinomeXplorer (Horn et al, 2014) was used for predicted kinase-substrate relations in phosphoproteomics clusters. TransFac (Matys et al, 2006;Wingender et al, 1996) database was used for transcription factor-target gene relations in transcriptomics clusters. Enrichment analysis of biological processes was done using the hallmarks gene sets from MSigDB (Liberzon et al, 2015).

Drug off assays
For long-term assay, WiDr cells were seeded in one 6-well plate (200,000 cells/well), starved in serum-free media for 24 h, and then treated with 3 µM PLX4032 and 3 µM gefitinib in complete medium. Treatment was interrupted after 78 days for a "drug off" period of 10 days. Pictures were acquired using an ECLIPSE Ti-e inverted microscope (Nikon, Tokyo, Japan) at a magnification of 10x. For short-term assay, WiDr cells were seeded in one 96-well plate (5000 cells/well) and 24 h after treated with 3 µM PLX4032 and 3 µM gefitinib. Treatment was interrupted after 5 days for a "drug off" period of 5 days. Cell growth was monitored in the IncuCyte™ automated microscope (Essen Bioscience, Ann Arbor, USA) and phase-contrast images were collected every 2 h using objective Nikon 10x. For complete information see Supplementary Methods.

Data and code availability
The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE (Vizcaíno et al, 2016) partner repository with the dataset identifier PXD007740. Sequence data has been deposited at the European Genome-phenome Archive (EGA) (Lappalainen et al, 2015)   In the phosphoproteomic analysis, selection of technical replicates involved discarding the poorest replicates before quantification. The resulting processed data were further filtered to only include data that was quantified in at least one timepoint of any condition. Analysis indicates good correlation (R>0.8) among the three biological replicates for all three omics datasets. C. Cross-omics Graphical User Interface output. The GUI enables exploration of the multi-omics data for specific genes under all the tested conditions. Users can select a measured gene of interest from a dropdown menu, and visualize complete data at the protein, phosphoprotein and transcript level, over the time-course and over all experimental conditions.
An example output for the gene CDK1 is shown. The user can toggle which number of sites is seen, but not the exact phoshposite.

Proteomics analyses
Lysed cells were defrost and each tube was supplemented with 300 µL of fresh lysis buffer.
Cells were further lysed by 10 rapid passages through 23G needle and by sonication on ice.
Cell debris were removed by centrifugation at 20,000 x g for 30 min at 4 °C and cleared supernatants were stored at -80 °C. The total protein concentration was measured using Bradford assay (Bio-Rad).
washed 3 times during 10 min with TBS-T (0.1 %). HRP-coniugated secondary antibodies (Bio-Rad) were diluted 1:10,000 in 5% BSA in TBS-T and incubated for 1 h at room temperature while shaking. Subsequently, membranes were washed additional 3 times during 10 min with TBS-T (0.1 %). Final protein detection was performed using Clarity ECL Western Blotting substrate (Bio-Rad) and blot imaging was performed using the Chemidoc Touch Imaging System (Bio-Rad).

Phosphopeptide enrichment
Ti 4+ -IMAC material was prepared and used as previously described (Zhou et al, 2013). Briefly, the affinity material was loaded onto GELoader tips (Eppendorf) using a C 8 plug. specified as enzyme and up to two missed cleavages were allowed. Cysteine carbamidomethylation was set as a fixed modification, while methionine oxidation and protein Nterm acetylation were set as variable modifications. Phosphorylation on serine, threonine and tyrosine was also selected as variable modification for the phosphoproteomics analysis. The allowed fragment mass deviation was set to 20 ppm for FTMS, and both minimum peptide length and the maximum peptide charge were set to 7. Fast Label free quantification (LFQ) was performed and 'match between runs' was enabled. Peptide and protein identification was set to 1 % FDR. The quantified output (proteinGroups.txt, phospho (STY)Sites.txt) were processed using a custom Python package (PaDuA). Potential contaminants and reverse peptides were removed, and filters for Class I phosphosites (localization probability > 75 %) and 'only identified by site' were applied for phosphoand protein data respectively. Normalization was performed by subtracting the median of log 2 transformed intensities from each column. For phospho-data, 'expand side table' function was applied before normalization. Median of technical replicates was performed for each dataset and resulting values were filtered to ensure each protein or phosphosite had valid measurements in at least one time point of any of the six cell culture conditions. Enrichment analysis was calculated using modificationSpecificPeptides.txt table.
Final processed output were exported for subsequent analysis in R.
containing drugs. In three 96 wells plates, gefitinib, lapatinib or sapitinib were serially diluted to a final concentration range of 120 nM to 30 µM as a single treatment and in combination with PLX4032 at a fixed concentration of 3 µM (n=4). In the fourth 96-well plate, PLX4032 was serially diluted to a final concentration range of 30 nM to 30 µM in combination with gefitinib at a fixed concentration of 3 µM and with either 1 mM DCA (n=4) or 10 µM etomoxir (n=4) respectively. In the last 96-well plate, PLX4032 was serially diluted to a final concentration range of 30 nM to 30 µM in combination with gefitinib at a fixed concentration of 3 µM (n=8).
After 96 h of treatment, 1 mM DCA was added to four replicates and 10 µM etomoxir was added to other four replicates. All plates contained a column with untreated cells as a reference sample and were treated until 7 days (37 °C, 5 % CO 2 ). Media containing the drugs was replaced after 72 h.
IC50 of single treatments was determined seeding WiDr cells in three 96-well plates at a density of 10,000 cells/well. After 24 hours of incubation (37 °C, 5 % CO 2 ), media was removed from all plates and replaced with media containing either PLX4032 or etomoxir or DCA serially diluted in four replicates to a final concentration range 30 nM-30 µM, 200 µM-0.2 µM and 100 mM-0.1 mM, respectively. All plates contained a column with untreated cells as a reference sample and were treated until 72 h (37 °C, 5 % CO 2 ).
Cell growth inhibition was monitored in all assays using the IncuCyte™ automated microscope (Essen Bioscience, Ann Arbor, USA) and phase-contrast images were collected every 2 h using a 10x Nikon objective. Phase confluence percentage from each well at each time point was exported into GraphPad Prism 7.0 software. The area under the curve (AUC) was calculated for each concentration (n=4), normalized in respect to untreated cells and fitted using a fourparameter logistic curve. Percentage of growth in single and combination treatments were visualized as dose response curves.
For long-term assay, WiDr cells were seeded in one 6-well plate (200,000 cells/well) and grown to around 60 % confluence. After 24 h starvation in serum-free media, cells were treated with PLX4032 and gefitinib both at a fixed concentration of 3 µM, in complete medium. Media was replaced two times per week and treatment was interrupted after 78 days. Pictures were acquired using an ECLIPSE Ti-e inverted microscope (Nikon, Tokyo, Japan) at a magnification of 10x. Cell confluence was measured using Fiji plugin in ImageJ software.
For short-term assay, WiDr cells were seeded in one 96-well plate at a density of 5000 cells/well. After 24 hours of incubation (37 °C, 5 % CO2), media was removed from all plates and replaced with media containing PLX4032 and gefitinib at fixed concentration of 3 µM. Media was replaced (once/twice) per week and treatment was interrupted after 5 days. Cell growth was monitored in the IncuCyte™ automated microscope (Essen Bioscience, Ann Arbor, USA) and phase-contrast images were collected every 2 h using objective Nikon 10x.