Introduction

PPINs are a major principle of biological organization. Several large-scale studies have demonstrated their importance in organizing fundamental cellular processes1,2,3,4. Increasing evidence suggests that PPINs are altered in human disease and that such changes contribute to pathogenesis3,5,6. However, we lack a clear comprehension of how such changes occur and how they affect PPIN function. This particularly applies to oncogenic mutations, where we rarely understand how system wide effects are generated. For instance, oncogenic RAS mutations, which activate the ability of RAS proteins to engage downstream effectors, occur in ~30% of all human cancers7. RAS mutated cancers are resistant to most targeted therapies, and inhibiting downstream effectors has proven ineffective, likely because of complex feedback structures in the downstream pathways and the large number of effector pathways8. These impediments highlight the need for a systems level understanding of the RAS signaling network and the changes associated with RAS transformation8. The KRASG13D mutation (mtKRASG13D) investigated here is the second most frequent RAS mutation in CRC7, and is associated with aggressive behavior and poor clinical outcomes9.

We use quantitative mass spectrometry (qMS) to map KRAS regulated PPINs in two closely related CRC cell lines that express either transforming or non-transforming levels of mtKRASG13D. Focusing on the epidermal growth factor receptor (EGFR) signaling network, where KRAS plays a key role in CRC10, we analyze 1710 immunoprecipitates and map >6000 PPIs involved in EGFR signaling (Fig. 1). To analyze this dataset we develop an analysis pipeline for the quantitative comparison of PPIN data between different cell lines. Finding that the expression of transforming levels of mtKRAS correlates with substantial rewiring of the EGFR PPIN, we analyze the functional consequences of this rewiring on protein complex assemblies, information flow, and biological responses including the prognosis of CRC patients. To facilitate the utilization of this extensive data for further research we develop PRIMESDB.eu (https://primesdb.eu/), an integrated database and analysis platform for exploring the PPINs described here.

Fig. 1: Experimental and data analysis workflow for the comparative mapping of PPIs in the EGFR network.
figure 1

Baits were chosen based on the core EGFR network described by Kiel et al.18 and additional manually curated literature information. Flag-tagged expression vectors were constructed using the Gateway cloning system and transfected into HCT116 (mtKRASHi) and HKE3 (mtKRASLo) cells. Careful titration of the transfected plasmids ensured similar protein expression in both cell lines grown in SILAC media as monitored by Western blotting. For MS experiments similar amounts of baits were expressed in SILAC labeled HCT116 and HKE3 cells and immunoprecipitated (IP) with anti-Flag antibodies. To assure robust quantitation the SILAC label was swapped, i.e. each bait was isolated from HCT116 and HKE3 cells grown in heavy or light medium, respectively. After trypsin digestion peptides were identified and quantified by orbitrap mass spectrometry. Raw data were analysed using MaxQuant59 and further processed using HiQuant19 implementing a stringent pipeline to retain only true interactors. Based on these data two quantitative protein–protein interaction networks, termed EGFRNetmtKRAS-Hi and EGFRNetmtKRAS-Lo, were reconstructed and are shown in a combined differential network representation. EV Ctrl, empty vector control transfection.

Full size image

Results

Mapping the EGFR PPIN in mtKRAS Cells

We mapped the effects of mtKRAS on PPINs downstream of the EGFR in HCT116 cells, which have been widely used to study mtKRAS functions in CRC. HCT116 harbor an oncogenic KRASG13D allele, which was targeted for disruption by homologous recombination to generate the non-tumorigenic HKE3 cells11. A thorough genetic, biochemical and biological characterization described in a previous publication12 and Supplementary Fig. 1 confirmed that HCT116 and HKE3 are closely related cell lines. Compared to HCT116 cells, HKE3 had a non-transformed phenotype, as reflected by EGFR inhibitor sensitivity and reduced migration, proliferation, colony forming ability, and anchorage independent growth (Supplementary Fig. 1). Interestingly, despite this non-transformed phenotype HKE3 cells retain a genetically stable low-level expression of mtKRASG13D, likely due to a duplication of the mutant mtKRASG13D allele in HCT116 and knockout of only one copy in the HKE3 cells (Supplementary Fig. 1A). Recent findings indicate that oncogenic KRAS mutations occur in normal tissues and that KRAS activity needs to exceed a threshold to drive cancer progression and metastasis13,14,15,16. Thus, this cell line pair offers the opportunity to compare the EGFR PPIN in cells expressing a transforming vs. a non-transforming dosage of mtKRAS. Furthermore, HCT116 are not addicted to mtKRAS for survival, minimizing selection pressure to acquire compensatory mutations when mtKRAS dosage is reduced17.

To attain a representative coverage of the EGFR PPIN we selected 95 bait proteins (Supplementary Data 1) based on a highly curated EGFR signaling network map18 and a literature survey of the EGFR pathways involved in CRC pathogenesis and progression. The baits cover the main functions of EGFR signaling including key kinases, phosphatases, scaffold and adapter proteins in the network. The baits were expressed as FLAG-tagged proteins carefully titrating transfection to achieve a similar expression level in both cell lines. Baits were immunoprecipitated and associated prey proteins were identified by high-resolution Orbitrap qMS. To ensure high data quality, we analyzed 95 bait and empty vector control immunoprecipitates (IPs) from both forward and reverse SILAC labeled cells using three biological and two technical replicates per bait resulting in 1710 samples and 1140 qMS analyses (Fig. 1). As common analysis methods for AP-MS data are ill suited for quantitatively comparing PPI data from different cell lines and from a biased bait selection, we used HiQuant19, which we specifically developed analyse MS data from complex experiments such as ours. The pipeline includes rigorous steps for data quality control, normalization, statistical analysis, and network construction. This workflow includes a stringent two-step procedure to exclude false positive interactors (Supplementary Methods, section 12). The EGFR PPINs in the HCT116 and HKE3 cells, termed EGFRNetmtKRAS-Hi and EGFRNetmtKRAS-Lo, respectively, were reconstructed from high-confidence bait–prey interactions.

EGFR PPIN network architecture

EGFRNetmtKRAS-Hi and EGFRNetmtKRAS-Lo consist of 3162 and 2788 bait-prey interactions, respectively (Supplementary Fig. 2A–E, Supplementary Data 2). This network size is within the expected distribution of known PPINs (Supplementary Fig. 3A). 93 of the 95 baits had a least one prey detected in both cell lines. More than 70% of preys in EGFRNetmtKRAS-Lo were also nodes in EGFRNetmtKRAS-Hi, indicating that most nodes are true components of the EGFR PPIN, since they were independently detected in both cell lines. Both EGFRNets are small-world, single-component networks (i.e., all nodes are reachable from one another) with similar scale-free topologies and comparable other network properties such as average path length, node degree, betweenness centrality (bc), and clustering coefficients (Supplementary Fig. 2F–I). For example, many of the major hubs (highly connected nodes) in EGFRNetmtKRAS-Hi are also hubs in EGFRNetmtKRAS-Lo including GRB2, RAB5A, RAF1, and SH2D3C. While GRB2 is a well-known hub coordinating different aspects of EGFR signaling20, some of the other high bc nodes are thought to have specialized functions. Our data suggest that these proteins are involved in a far greater degree of crosstalk with other cellular processes than previously assumed.

A comparison with currently known PPIs shows that >80% of interactions discovered in our study are new, attesting to the value of focused PPIN mapping studies complementing genome-wide efforts and suggesting that many PPIs may be highly dependent on the cellular context. Such a high proportion of novel interactions is consistent with other large-scale, AP-MS based, interactome mapping efforts1,2,4,21,22. Testing 17 arbitrarily chosen interactions using conventional co-immunoprecipitation/Western blot experiments showed that the PPI data were highly reproducible by a different method (Supplementary Fig. 3B, C).

As AP-MS may underrepresent interactors of integral membrane proteins23, we used MYTH, a membrane yeast two-hybrid assay24, to identify binary protein interactors of the human EGFR family, ERBB1–4 (Supplementary Methods, section 30). This interactome map comprised 405 interactions (Supplementary Fig. 4 and Supplementary Data 3) including 181 new interactions. All bait–prey interactions detected in theEGFRNets can be explored in the Supplementary Data and at http://primesdb.eu/.

mtKRASHi induces extensive PPIN rewiring

To identify interactions that were significantly rewired in mtKRASHi cells, i.e., interactions that were gained/lost in one EGFRNet or present in both but with significantly altered prey abundance, we statistically compared prey abundance between each bait-prey complex in the EGFRNets. Of the 4420 bait–prey interactions detected in at least one EGFRnet (Supplementary Data 4), 1368 were significantly rewired i.e., prey abundance was significantly different between the two PPINs at P ≤ 0.05, significance A ≤ 0.05 (Supplementary Data 5). Six hundred and thirty four of the rewired interactions were edges only in EGFRNetmtKRAS-Hi, and 406 were edges only in EGFRNetmtKRAS-Lo indicating that most rewiring is due to interaction gains or losses (Fig. 2). The 328 remaining rewired interactions were present in both networks but with significantly different prey abundances (P ≤ 0.05, significance A ≤ 0.05). These data suggest that the oncogenic mtKRAS activity in mtKRASHi cells initiates a ripple effect throughout the network substantially altering network topology far beyond direct KRAS interactors.

Fig. 2: The EGFRNetmtKRAS-Hi and EGFRNetmtKRAS-Lo PPINs are rewired.
figure 2

a The number of preys identified for each bait-prey AP-MS complex. Red, rewired preys enhanced in mtKRASHi cells; blue, rewired prey proteins enhanced in mtKRASLo cells. AP-MS complexes are named based on the bait protein and are shown on the x-axis. b Network spoke model view of rewired interactions. Bait–prey interactions identified in both networks with similar prey abundance are shown in gray. Bait–prey interactions that were identified only in EGFRNetmtKRAS-Hi or EGFRNetmtKRAS-Lo are shown as solid red or blue edges, respectively. Bait–prey interactions that were detected in both networks but where prey abundance was significantly higher in EGFRNetmtKRAS-Hi or EGFRNetmtKRAS-Lo are shown as dotted red or blue edges, respectively. Visit primesdb.eu to explore the complexes in more detail/with greater resolution. c Zoom-in on four nodes that represent mixed, preferential or no rewiring. Source data are provided as a Source Data file.

Full size image

Potential drivers of PPIN rewiring

To investigate which molecular mechanisms could drive PPIN rewiring, we first analyzed whether genetic mutations other than KRASG13D played a role, since genetic variation has previously been associated with PPIN rewiring25. Using whole genome sequencing we identified genetic alterations, including copy number variations (CNVs), insertions/deletions (InDels), synonymous and nonsynonymous single-nucleotide-variants (SNVs) between the two cell lines (Supplementary Data 6–8; Supplementary Fig. 5A). Using the Genome Analysis Toolkit26 27 genes were predicted to be impacted by structural variants, but no gene was a node in the EGFRNets. Considering CNVs, five genes were EGFRNet nodes, but only one gene product, PPP3CA, was rewired. Of the 170,135 SNVs and small InDels found different between mtKRASHi and mtKRASLo cells 1091 were variants of predicted high/medium impact27 (Supplementary Data 6). Of these, 70 were nodes in the EGFR PPI network and 36 were rewired. Considering that EGFRnets contain 4420 nodes, of which 1360 have rewired interactions, SNVs affect 1.6% of nodes and 2.6% of rewired interactions. These data suggest that structural variants, SNVs and CNV-driven changes in gene/protein expression have limited impact on EGFRNet rewiring. Nonetheless, we cannot rule out that these or other genetic differences influence some PPIs by affecting gene promoter usage, mRNA editing, or codon usage. We also considered that rewired prey could simply represent lowly or highly expressed nodes. However, we found no bias in the gene expression distribution of rewired nodes compared to unchanged nodes (Supplementary Fig. 5B) suggesting that genetic changes that alter gene/protein expression, e.g., CNVs, do not make major contributions to PPIN rewiring.

To further explore this, we directly tested whether changes in protein expression between the two cell lines are linked to the observed EGFRNet rewiring. We profiled protein abundances in the mtKRASHi and mtKRASLo cell lines using qMS (Supplementary Data 9). 404 of the 4685 proteins quantified showed a significant difference in abundance (P ≤ 0.05). Pathway analysis revealed that proteins more abundant in mtKRASHi were enriched for roles in the cell cycle, consistent with the increased proliferation rate of these cells (Supplementary Fig. 1E). By contrast, proteins more abundant in mtKRASLo cells were enriched for roles in oxidative phosphorylation, lipid metabolism, and the lysosome. The decreased expression of proteins involved in oxidative phosphorylation in mtKRASHi cells is consistent with a metabolic switch from oxidative phosphorylation to glycolysis, a hallmark of cancer cells known as the Warburg effect28. A strong relationship between KRASG12D dosage and increased glycolysis was recently reported29. Similarly, lipid metabolism reprogramming is also a hallmark of cancer cells, including CRC cells30.

We found a weak (r2 = 0.18) but significant correlation (P < 0.001) between fold-change in abundance in the AP-MS protein complexes and fold-change in protein expression between the cell lines (Supplementary Fig. 5C). One hundred and fourteen differentially expressed (DE) proteins were nodes in the EGFRNets (Fig. 3a, Supplementary Fig 5D, E), two of them corresponding to baits (RPS6KA1 and SH3KBP1, Δ = 1.7-fold). This was not more than statistically expected (P = 0.054) indicating that the EGFR network was not especially enriched for DE proteins. However, 74 of the 114 DE proteins represented rewired nodes, which was statistically significant (P = 4.21E−5), confirming an association between differential node abundance and network rewiring. Interestingly, some bait-prey complexes were particularly enriched for DE proteins (Supplementary Fig. 5F, G). For example, of 71 rewired preys in the SH2D3C complex, 16 (22%) were also DE. Overall, these data suggest that differences in protein expression between mtKRASHi and mtKRASLo cells may underlie some of the rewired interactions. However, this association was lost when considering proteins at DE >2-fold. Furthermore, as ~90% of rewired nodes were not identified as DE proteins, differences in protein expression alone cannot explain the wide-spread network rewiring. We have not investigated the reverse possibility that PPIs may affect protein stability31,32.

Fig. 3: Potential drivers of the EGFR PPI network rewiring.
figure 3

a The number of rewired prey proteins for each bait-prey AP-MS complex that were assessed for differential protein expression between the two cell lines. Rewired prey proteins that were significantly more abundant in the mtKRASHi cells are shown in red. Rewired prey proteins that were significantly more abundant in the mtKRASLo cells are shown in blue. Four selected AP-MS complexes highlighting differentially abundant prey proteins (larger nodes) are also shown. b The number of rewired prey proteins for each bait-prey AP-MS complex that were assessed for differential phosphorylation between the two cell lines. Rewired prey proteins that were significantly more phosphorylated in the mtKRASHi cells are shown in red. Rewired prey proteins that were significantly more phosphorylated in the mtKRASLo cells are shown in blue. Four selected AP-MS complexes highlighting differentially phosphorylated prey proteins (larger nodes) are also shown. c Statistically enriched pathways among the 735 prey proteins involved in rewired interactions. Source data are provided as a Source Data file.

Full size image

To assess other potential drivers of PPIN rewiring, we examined protein phosphorylation, which can generate docking sites and affect binding affinities between proteins, thereby influencing protein complex formation. qMS-based phosphoproteome analysis of mtKRASHi and mtKRASLo cells identified 384 differentially phosphorylated proteins (Supplementary Data 10). Two hundred and seventy one proteins were preferentially phosphorylated in mtKRASHi cells and were enriched for roles in cell cycle and apoptosis related pathways, while the 121 proteins preferentially phosphorylated in mtKRASLo cells were weakly enriched for cytokine signaling and related processes (Supplementary Data 10). Eighty nine differentially phosphorylated (DP) proteins were nodes in the EGFRNets (Fig. 3b). Compared to the number of network nodes that were also represented in the phosphoproteomics screen, this was not more than statistically expected (P = 0.06) indicating that the EGFRNets were not enriched for DP proteins. However, 56 of the 89 DP proteins (63%) mapping to EGFRNets were also significantly rewired nodes (P < 0.01). Interestingly, this association was even stronger when considering interactions that were enhanced (or only found) in EGFRNetmtKRAS-Hi. Rewiring of these interactions was significantly correlated with higher phosphorylation in the mtKRASHi cells (P < 0.001). These results suggest that differential phosphorylation could contribute to PPIN rewiring, especially when mtKRAS signaling is high (Supplementary Fig. 5F, G). However, as with differentially expressed proteins, the majority of rewired nodes were not associated with differential phosphorylation. This suggests that the observed network rewiring is an emergent property of the changed cellular state in cells expressing transforming levels of mtKRAS and is not readily predicted by changes in any single factor such as protein expression or phosphorylation.

PPIN rewiring modifies protein complexes and their functions

Pathway analysis of all 735 prey proteins involved in rewired interactions revealed a statistically significant enrichment for roles in processes including RNA splicing, mitochondrial translational, protein folding by the chaperonin-containing-TCP1 (CCT) complex and cell migration (Fig. 3c, Supplementary Data 11). Interestingly, CRISPR and shRNA-based screens of HCT116 cells and isogenic wild-type KRAS derivatives found that synthetic lethal KRASG13D genes had roles in mRNA splicing and mitochondrial translation, and that these processes were required for KRASG13D oncogenicity33. Proteins encoded by synthetic lethal genes identified in this CRISPR screen were significantly enriched (9 of 55; P < 0.01) for rewired PPIs in our network. The shRNA screen in this study identified several CTT components as synthetic lethal genes. In our study, several baits co-precipitated the entire CCT complex in both mtKRASHi and mtKRASLo cells (Supplementary Data 4), with interactions between CCT and the baits BMX, LCK, and PTK6, being significantly rewired (Supplementary Data 5).

Similar correlations were observed when assessing how PPIN rewiring affected the composition of known protein complexes described in CORUM, a curated database of experimentally determined mammalian protein complexes34. We detected 42 and 40 CORUM complexes where at least 70% of their component proteins were nodes in EGFRNetmtKRAS-Hi and EGFRNetmtKRAS-Lo, respectively (Supplementary Fig. 6, Supplementary Data 12). Several of these CORUM complexes are involved in processes mediated by EGFR signaling including actin cytoskeleton organization, RAF, and NFκB signaling35,36. However, other complexes participate in functions not usually associated with EGFR signaling, e.g. chromatin modification, regulation of protein folding, mRNA splicing and protein translation. This analysis suggests that PPIs organize different aspects of EGFR signaling including roles that have not been characterized yet. Most CORUM complexes were present in both networks, although some were extensively rewired. CORUM complexes, where >60% of constituent proteins were rewired preys, included complexes involved in mRNA transcription, splicing and protein folding (Supplementary Fig. 6) suggesting that such house-keeping functions support mtKRAS transformation. Furthermore, complex formation can stabilize proteins31 and may contribute to the differential protein abundance between mtKRASHi and mtKRASLo cells.

Interestingly, rewired interactions were non-randomly distributed across the bait-prey complexes (Supplementary Fig. 7). Many bait-prey complexes, including AKT1, FGR, HDAC1, LYN, MAPK3 (ERK1), RIPK1, and TNK1 complexes, had no significantly rewired interactions, while others, such as MAP2K1 (MEK1), RAC1, and SH2D3C, were substantially rewired (Supplementary Data 5). In some cases, rewiring predominantly involved the gain of new bait–prey interactions in mtKRASHi cells. For example, all 25 rewired bait–prey interactions in the BAD complex were only detected or significantly more abundant in the EGFRNetmtKRAS-Hi (Fig. 4a). BAD is a proapoptotic protein37, which contributes to the higher apoptosis rate of mtKRASHi vs. mtKRASLo cells, as shown by siRNA knockdown experiments (Fig. 4b). Interestingly, the BAD interactions enhanced in EGFRNetmtKRAS-Hi included a higher abundance of protein phosphatase 2A (PP2A) family members (Supplementary Data 5), which correlated with a lower phosphorylation of BAD on S112 and S155 in HCT116 (Fig. 4c). These sites inactivate BAD’s pro-apoptotic function and can be phosphorylated by cAMP dependent protein kinase (PRKA)37, which interacted with BAD in both cell lines. Inhibition of PRKA reduced S112 and S155 phosphorylation in mtKRASLo cells, while PRKA inhibition increased BAD phosphorylation in mtKRASHi, especially on S112 (Fig. 4c). This differential action of PRKA is likely due to its ability to activate PP2A by phosphorylating the B56δ subunit (PPP2R5)38, which is mainly bound to BAD in mtKRASHi cells. Consequently, PRKA inhibition preferentially enhanced apoptosis in mtKRASLo cells (Fig. 4d). These results suggest that PPI rewiring can profoundly subvert the biological effects of PRKA signaling, in this case converting a survival signal into a pro-apoptotic signal (Fig. 4e).

Fig. 4: Differential regulation of BAD protein phosphorylation and its biological effect.
figure 4

a PPI interactions in the BAD complex. Red broken lines, PPIs enhanced in EGFRNetmtKRAS-Hi vs. EGFRNetmtKRAS-Lo; red solid lines, EGFRNetmtKRAS-Hi exclusive interactions; gray, unchanged interactions. b Knocking down BAD expression significantly reduces apoptosis in mtKRASHi but not mtKRASLo cells. Apoptosis was measured 24 h post treatment. Ctrl., untargeted siRNA; BADkd, BAD specific siRNA. The reduction of BAD protein expression was assayed by Western blotting using tubulin (Tub) as loading control. Apoptosis assays represent three independent experiments; error bars are SD, and * means significant (P < 0.05) according to Student’s t-test; ns, not significant. c mtKRASHi and mtKRASLo cells were treated with the PRKA inhibitor H89 as indicated before BAD phosphorylation at S112 and S155 were assessed by Western blotting using phospho-specific antibodies. Numbers below lanes represent BAD phosphorylation normalized to total BAD protein expression. Samples shown in b and c are from the same Western blots, where irrelevant lanes were removed as indicated by vertical lines. d Apoptosis in response to 20 μM H89 treatment was assessed as in a. The data represent two independent experiments with 3 and 4 biological replicates, respectively. Error bars are SD, and * means significant (P < 0.05) according to Student’s t-test; ns, not significant. e A proposed model of the differential biological effect of PRKA due to differential PPIs. Source data are provided as a Source Data file.

Full size image

Another substantially rewired node was PTK6 (Supplementary Fig. 8). PTK6 is a poorly characterized tyrosine kinase, which is amplified or overexpressed in 16% of CRC patients. PTK6 can stimulate CRC cell survival and oncogenic signaling in a kinase dependent manner, but suppresses epithelial-to-mesenchymal transition in a kinase independent fashion39. PTK6 rewiring mostly decreased interactions with the CCT chaperonin complex in mtKRASHi cells, whereas the interaction with metastasis associated 1 family member 2 (MTA2) was increased (Supplementary Fig. 8A). Given the key role that MTA2 plays in cell motility and metastasis40, we examined whether PTK6 contributed to the differential cell migration observed between mtKRASHi and mtKRASLo cells (Supplementary Fig. 1F). mtKRASHi and mtKRASLo expressed similar amounts of endogenous PTK6 (Supplementary Fig. 8B). Overexpressing PTK6 accelerated migration in HCT116 cells but inhibited it in HKE3 cells (Supplementary Fig. 8C). Increased migration was dependent on PTK6 kinase activity, as a kinase dead PTK6 mutant slowed migration (Supplementary Fig. 8D). Knocking down endogenous PTK6 decreased migration specifically in mtKRASHi but not in mtKRASLo cells (Supplementary Fig. 8E). Taken together, these results suggest that PTK6 preferentially enhances migration in cells with high KRAS activity.

Network rewiring alters information flow through EGFRNets

The extensive changes in network wiring and protein complex composition suggested that PPIN rewiring may alter signal processing in the EGFR network. First, we explored how different concentrations of active KRAS affects the formation of KRAS complexes with known effector proteins. Activated RAS proteins signal by binding a range of effectors through a single, shared binding domain7 leading to competition between effectors. To analyze the formation of specific KRAS-effector complexes we constructed an equilibrium binding model of proteins competing for a single target (see Methods section). This model classifies KRAS effectors into low and high affinity binders, whose binding dissociation constants (Kd’s) are greater or smaller, respectively, than the abundance of active KRAS. It shows that for low-affinity effectors the corresponding KRAS complex concentrations are proportional to the effector concentration divided by the Kd, whereas for high-affinity interactors the resulting KRAS complexes concentrations are determined by the abundance of active KRAS and effectors alone. Thus, changes in mtKRAS concentration can profoundly rearrange the composition of KRAS-effector complexes, which rather than changing the strengths of downstream pathway activation shifts signaling from high to low affinity effectors as mtKRAS dosage increases (Fig. 5a). As measured by quantitative Western blotting, the mtKRAS concentrations in mtKRASLo and mtKRASHi cells are ~150 nM and ~400 nM, respectively, indicating that high affinity RAS effector complexes prevail in mtKRASLo cells, while low affinity effectors dominate signaling in mtKRASHi cells. Specifically, the model predicted that fold-changes in KRAS-bound fractions are higher for low-affinity than for high-affinity effectors. Ranking effectors by the fold-changes in KRAS-bound fractions allowed us to estimate their relative contribution to downstream signaling. Applying this analysis to baits that participate in bona-fide KRAS effector pathways (RAF/MAPK, RAL, PI3K, TIAM, AFDN, PLCε, and RIN1), we calculated the sensitivity of a node responding to different mtKRAS doses by summing the log fold-changes of interactions (normalized by the number of pathway nodes measured) in each bait-prey AP-MS complex. These experimentally deduced sensitivity ranks of KRAS-effector complexes correlated with the model-predicted ranks (Supplementary Data 13). These results suggest that the threefold difference in mtKRAS activity induces the formation of very different KRAS-effector complexes that initiate network rewiring by engaging different signaling pathways (rather than stronger activate the same set) that propagate changes further downstream.

Fig. 5: KRAS effector pathway and information flow (IF) analysis of the EGFRNetmtKRAS-Hi and EGFRNetmtKRAS-Lo networks.
figure 5

a Dependence of KRAS-effector complex concentrations of effectors binding with high (RAF1, blue) or low affinity (RALGDS, red) on the abundance of mtKRAS. Broken gray lines indicate KRAS concentrations in mtKRASLo and mtKRASHi cells. b Plot showing nodes in the top 5th percentile in terms of their IFS and predicted to receive more flow in EGFRNetmtKRAS-Hi (red) or EGFRNetmtKRAS-Lo (blue). c Transcription factors (TFs) with at least 20% higher information flow in EGFRNetmtKRAS-Hi (red) and EGFRNetmtKRAS-Lo (blue). Gene expression of d FOXO1, e MYC, f FOS, and g STAT1 as determined by RNAseq analysis of TGFα stimulated cells. ***EdgeR FDR < 0.001; ****EdgeR FDR < 0.0001. The top five enriched transcription factor binding (TFBS) site motifs in the promoters of genes upregulated in h mtKRASHi and i mtKRASLo cells. j Reporter gene assays of the activity of STAT1/2 and STAT3. Error bars represent standard deviation, and P values in K were determined by a two-tailed Student’s t-test. *P < 0.05; **P < 0.01. The data represent three independent experiments. Source data are provided as a Source Data file. Mathematica code for 5A is provided in Supplementary Software 1.

Full size image

Given these extensive changes in network wiring and protein complex composition we hypothesized that PPIN rewiring would also alter signal processing, leading to differential activation of downstream transcriptional programs. In order to investigate this hypothesis in an unbiased way not limited to known KRAS effector pathways, we employed a computational modeling based approach called information flow (IF) analysis41,42. This method simulates IF in a network through discrete-time random walks from a source node, i.e., EGFR, to downstream sinks, i.e., transcription factors (TFs). To model the impact of PPIN rewiring, we simulated IF independently in the EGFRNetmtKRAS-Hi and EGFRNetmtKRAS-Lo networks and calculated an IF score (IFS) for each node in the two networks that reflects the volume of signals flowing through a node. Nodes with high IFS in both networks included known key transducers of EGFR signaling, e.g., GRB243, indicating that major hubs are used regardless of mtKRAS dosage. However, 119 nodes had a >2-fold difference in IFS in EGFRNetmtKRAS-Hi vs. EGFRNetmtKRAS-Lo (Fig. 5b and Supplementary Data 14), indicating potentially critical differences in signal processing. Interestingly, many of the highest scoring nodes that received more information flow in the EGFRNetmtKRAS-Hi network were proteins involved in protein folding including heat shock protein (HSP) 70 family members (HSPA1A and HSPA8), and HSP90 family members (HSP90AA1, HSP90AB1, and HSP90AB3P). HSP70 and HSP90 expression is upregulated in many cancers including CRC44, and high HSP70 expression is associated with poor clinical outcomes in CRC45. Furthermore, HSP90 inhibitors are in clinical trials for several cancers including CRC46. Another node with higher IFS in EGFRNetmtKRAS-Hi was SRC, which is a major promoter of CRC proliferation, metastasis, drug resistance, and is overexpressed in ~80% of CRCs47. These data suggest that network nodes with increased IF in mtKRASHi cells contribute to the molecular pathogenesis of CRC and may represent potential drug targets.

Next, we assessed whether PPIN rewiring alters IF to transcription factors (TFs) in the EGFRNets (Fig. 5c). FOXO1 and MYC were predicted to receive higher IF in EGFRNetmtKRAS-Hi. Assessing gene expression in both cell lines by RNAseq prior to and following EGFR activation by TGFα, revealed that FOXO1 and MYC were more highly expressed in HCT116 cells (Fig. 5d, e). On the other hand, TFs including STAT1 and FOS received higher IF in EGFRNetmtKRAS-Lo, and their gene expression was significantly elevated in HKE3 cells (Fig. 5f, g). These results were consistent with the IF model predictions. Next, we analyzed TF binding sites in the promoters of genes that were differentially regulated between mtKRASHi and mtKRASLo cells. Genes upregulated in mtKRASHi cells were enriched for MYC binding sites (Fig. 5h), consistent with the prediction that MYC receives more IF through EGFRNetmtKRAS-Hi. Conversely, genes upregulated in mtKRASLo were enriched for the interferon-stimulated response element (ISRE) motif (Fig. 5i), a key motif in the promoters of STAT1/2-regulated genes48. The difference in FOS gene expression between the cell lines was particularly evident 60 min post-TGFα stimulation. Consistent with the prediction of higher FOS regulation through EGFRNetmtKRAS-Lo, the AP-1 binding site motif for FOS/JUN dimers was enriched in the promoters of genes upregulated in mtKRASLo cells at this timepoint (data not shown). The prediction that STAT1 receives lower IF in EGFRNetmtKRAS-Hi is consistent with reports that mtKRAS inhibits STAT1/2 expression49. To directly examine STAT activity, we used luciferase reporter genes that are regulated by STAT1/2/3 TFBS (Fig. 5j). STAT1/2 reporter activity was significantly elevated in mtKRASLo, while STAT3 activity was similar in mtKRASHi and mtKRASLo cells.

In summary, these data suggest that mtKRAS mediated PPIN rewiring alters signal flow through the EGFR network leading to the induction of different transcriptional programs. These analyses also support our PPIN reconstruction and the functional consequences of PPIN rewiring by unbiased global approaches.

Alterations in rewired baits predict CRC patient survival

The results presented above suggest that PPIN rewiring is associated with mtKRAS signaling and oncogenic potential. Therefore, we investigated whether alterations in bait proteins showing the most rewiring were prognostic of CRC patients’ clinical outcomes. We assessed survival data of 629 CRC patients from the TCGA dataset50. Fifty-four percent of patients had genetic or expression alterations in the top 20 most rewired bait proteins, as defined by the sum of rewired interactions associated with each bait (Fig. 6a, Supplementary Data 15). Patients with alterations in top 20 most rewired baits had significantly poorer survival (P < 0.04) than patients without alterations in these proteins (Fig. 6b). Ten-year survival was 34.61% for patients with alterations in the top 20 rewired baits vs. 61.43% for patients without. These data were robust to removing the bottom 50% of least significant (based on the significance A value) rewired interactions from the rewiring analysis and recalculating the top 20 most rewired baits. 18 of the 20 original top 20 baits were the same in this analysis (data not shown). By contrast, there was no significant association between alterations in the 20 least rewired bait proteins (Fig. 6c) and patient survival (P = 0.20) (Fig. 6d), although all of these baits were preselected because of their roles in the EGFR pathway. This association with survival became even stronger (P = 9.855e−3), if we defined the top 20 most rewired baits based on interactions that were selectively enhanced in the mtKRASHi cells (Supplementary Fig. 9A).

Fig. 6: PPIN rewiring and CRC prognosis.
figure 6

a The top 20 most rewired bait proteins. Interactions where the prey protein was identified only in EGFRNetmtKRAS-Hi or EGFRNetmtKRAS-Lo are shown as solid red or blue lines, respectively. Rewired bait-prey interactions where prey abundance was significantly higher in EGFRNetmtKRAS-Hi or EGFRNetmtKRAS-Lo are shown as dotted red or blue lines, respectively. Bait–prey interactions which were not significantly different are gray. b Six hundred and twenty-nine CRC patients from the TCGA were divided into two groups, those with alterations in the top 20 rewired bait proteins (339; 54%) and those without alterations in the top 20 rewired bait proteins (290; 46%). The alterations assessed were mutations, copy number changes, mRNA expression changes, and protein expression changes. Kaplan–Meier survival curves were plotted for the two patient groups using PRISM 7.0.3. Five-year survival was 53.5% for patients with genetic alterations affecting the top 20 rewired nodes compared to 68.5% for patients without alterations in these proteins, and ten-year survival was 34.61 vs. 61.43%, respectively. The log-rank test was used to assess statistical significance. c The bottom 20 least rewired bait proteins. d There was no significant (NS) difference in survival between patients with alterations in the bottom 20 least rewired bait proteins and those without. Source data are provided as a Source Data file.

Full size image

To assess the accuracy of the top-20 bait proteins to classify patients into high and low risk groups, we trained a Lasso classifier51 using the CRC patient data from TCGA. Five-fold cross-validation by subsampling the patient data into training (80%) and test (20%) datasets gave an accuracy of up to 0.79 (mean 0.70) and an area under the ROC curve (AUC) of 0.763 (Supplementary Fig. 9B). A similar classification using the bottom 20 least rewired proteins gave a much lower mean accuracy of 0.4 and AUC of 0.522. Several top rewired baits were highly connected nodes. Therefore, to ascertain that the association with patient outcomes was due to the rewiring of these baits and not just because they were highly connected, we selected the bottom 36 least rewired baits that together accounted for at least the same number of interactions as the top 20 baits. Patients with alterations in the top 20 rewired baits again showed significantly poorer survival (P < 0.017, log-rank test) than patients with alterations in these baits. (Supplementary Fig. 9C). Patients with alterations in the top 20 rewired baits also showed significantly poorer survival after adjusting for age and tumor stage (P < 0.03, log-rank test). These data suggest that the proteins, which we found to be the most rewired in mtKRASHi cells, are clinically relevant as alterations in these proteins are prognostic of CRC patient outcomes.

Discussion

Global PPIN mapping has validated the concept that the cell organizes its proteome as modules of PPIs that enable it to carry out its specific biological functions1,2,4. Many disease-associated mutations affect PPIs25, but the extent of adaption to disease mutations at a PPIN level and its functional consequences are unknown. Our comparative mapping of >6000 PPIs in the EGFR network in cells with low and high mtKRAS signaling reveals a widespread rewiring of the EGFR signaling network. Interestingly, rewiring percolates through the whole network and alters interactions that occur between core components of the EGFR pathway as well as interactions between proteins involved in downstream and seemingly peripheral processes. This suggests that enhanced mtKRAS activity results in extensive adaptive changes that are reflected by a reorganization of the PPIN. Genetic mutations associated with disease often alter PPIs25. However, the deep network propagation of PPI changes arising from a single mutation was unexpected and may explain why blocking mtKRAS signaling by inhibiting single effector pathways is ineffective8. Our global analysis and validation of the functional consequences of PPIN rewiring facilitates the rational design of combinatorial targeting of mtKRAS effectors, especially as PPIs gained in mtKRASHi cells inversely correlate with CRC patient survival. For instance, HSPs receive high IF in mtKRASHi cells, and HSP90 inhibitors recently were found to enhance the effects of conventional CRC drug therapies52. Likewise, our results that phosphorylation changes contribute to PPIN rewiring may “repurpose” kinase inhibitors as PPIN rewiring agents. We recently showed that mtKRAS also profoundly changes the metabolic and transcriptional landscapes of CRC cells53 confirming that mtKRAS widely affects cellular regulation on different levels. Surprisingly, our analysis shows that a low-level expression of mtKRAS is compatible with normal (i.e., untransformed) biochemical and biological behavior. This finding suggests that oncogenic mutations must reach a threshold activity before they produce a pathogenetic phenotype. While this view challenges the genetic mutational dogma of carcinogenesis, it reconciles with recent data finding oncogenic mutations in normal tissues54,55. The easy accessibility and analysis of our PPIN data through PRIMESDB.eu and DyNet56, an application for the visualization and analysis of dynamic molecular interaction networks will support systematic efforts of combinatorial mtKRAS pathway targeting.

While our study is comprehensive and integrates PPIN reconstruction with computational model analysis of network functions, it also has limitations. Although we have thoroughly characterized the two cell lines by WGS, biochemical and biological assays, we cannot formally exclude that differences other than mtKRAS activities contribute to our results, e.g., differential bait expression vs. endogenous levels, epigenetic differences, or nonsynonymous mutations that affect splicing or codon usage. The influence of such factors could be addressed by reconstitution experiments that titrate mtKRAS dosage and by studying other isogenic cell line pairs. Investigating all these aspects was beyond the scope of this study. However, we assessed and found no statistical association between rewiring and alternative splicing (data not shown), which in binary interaction screens substantially changed PPIs57. It also will be important to test whether the observed PPI rewiring is common to different mtKRAS cancer types. Our findings that PPI changes correlate with CRC patient prognosis and often affect proteins that are synthetic lethal with KRASG13D33 indicate that a core signature of consistently altered PPIs may exist in mtKRAS cells.

In summary, these results suggest that dynamic PPIN adaptations play major roles in translating the effects of genetic mutations into quantitative functional effects that re-direct information flow through signaling networks and reprogram biological outcomes.

Methods

Cell lines and cell culture

HCT116 (mtKRASHi) and HKE3 (mtKRASLo) cells11 were provided by Drs Shirasawa and Sasazuki. Cell lines were authenticated by whole genome sequencing (Supplementary Datas 6–8) and RNAseq as recently described12.

Baits and expression vectors

For AP-MS experiments 95 baits (Supplementary Data 1) were selected that provide a broad coverage of the known EGFR signaling network. Bait cDNAs were obtained from Origene and cloned into the SF-TAP vector58 with the FLAG-tag at the N-terminus using the Gateway cloning system (Thermo Fisher).

AP-MS experiments

Cells were transfected with baits titrated to achieve similar expression levels between the cell lines and labeled with stable isotope (SILAC) medium (Fig. 1). Baits were immunoprecipitated with anti-FLAG-M2 conjugated agarose beads (Sigma-Aldrich A2220), digested with trypsin and analyzed by quantitative MS using a Q-Exactive mass spectrometer (Thermo Fisher Scientific). Data were analyzed with MaxQuant59 and HiQuant19 software packages. See Supplementary Methods, sections 8–12, for a detailed description.

Protein expression profiling

Lysates of SILAC labeled cells were digested with trypsin and analyzed on an Orbitrap Fusion Tribrid mass spectrometer (Thermo Fisher Scientific). Data were analyzed with MaxQuant59 and HiQuant19 software packages. See Supplementary Methods, sections 13–15, for a detailed description.

Phosphoproteomics

Cell lysates were digested with trypsin, phosphopeptides were enriched using TiO2 beads and analyzed on a Q-Exactive mass spectrometer (Thermo Fisher Scientific). MS data were analyzed with MASCOT 2.4 (Matrix Science Ltd) and Progenesis 4 (Nonlinear Dynamics). See Supplementary Methods, sections 16–18, for a detailed description.

Construction of protein–protein interaction networks (PPIN)

The EGFRNetmtKRAS-Hi and EGFRNetmtKRAS-Lo networks were separately constructed by combining bait–prey interactions from each of the 95 chosen baits. Bait-prey interactions were included in the networks, if the abundance of the prey protein in the pull-down was significantly higher (P ≤ 0.05) than in empty vector controls and the significance A value for the prey protein was also ≤0.05. To identify interactions that were significantly “rewired” in the HCT116 (EGFRNetmtKRAS-Hi) network compared to the HKE3 (EGFRNetmtKRAS-Lo) network, we used HiQuant to directly compare the SILAC data from the two cell lines. We defined interactions as being “rewired” in EGFRNetmtKRAS-Hi, if the abundance of the prey protein in the pull-down was significantly different (P ≤ 0.05) compared to EGFRNetmtKRAS-Lo and the significance A value for the prey protein was also ≤0.05. Interactions that were identified in only one EGFRNet, but where prey abundance was subsequently not found to be statistically significantly different in the respective bait–prey complexes in the two cell lines were not considered as rewired interactions. The top rewired nodes in EGFRNetmtKRAS-Hi were identified as those with most rewired interactions. The topological properties of the networks including node degree, betweenness centrality, clustering coefficient and network scale-freeness were analyzed using the NetworkAnalyzer application60 in Cytoscape 361. Cytoscape session files for the EGFRNetmtKRAS-Hi, EGFRNetmtKRAS-Lo networks can be provided upon request.

Protein abundance and phosphorylation enrichment analysis

To investigate whether nodes in the EGFRNetmtKRAS-Hi network were enriched for differentially abundant proteins, a hypergeometric test was performed with the following parameters:

$$p\left( {X\, \ge k} \right) = \mathop {\sum }\limits_{x\,=\,k}^n \,\frac{{\left( {\begin{array}{*{20}{c}} K \\ x \end{array}} \right)\left( {\begin{array}{*{20}{c}} {N - K} \\ {n - x} \end{array}} \right)}}{{\left( {\begin{array}{*{20}{c}} N \\ n \end{array}} \right)}},$$

where N is the total number of proteins assayed in the protein expression analysis, n the total number of differentially abundant proteins identified, K the number of proteins in the EGFRNetmtKRAS-Hi network that were assayed in the protein expression analysis, k the number of differentially abundant proteins observed in the EGFRNetmtKRAS-Hi network.

A similar analysis was conducted to determine whether rewired nodes were enriched for differentially abundant or phosphorylated proteins. See Supplementary Methods, section 19, for a detailed description.

Equilibrium binding model of RAS binding partners to RAS-GTP

In order to determine how the concentrations of KRAS-effector complexes change with the concentration of active RAS in mtKRASHi and mtKRASLo cells, we developed a dynamic mathematical model that allowed us to investigate how competition for the single effector binding site on RAS and different abundances of low and high affinity RAS effectors impact the formation of RAS-effector complexes. See Supplementary Methods, sections 22–25, for a detailed description.

Information flow analysis of EGFRNetmtKRAS-Hi and EGFRNetmtKRAS-Lo networks

In order to analyze how the EGFR PPINs transduce information we employed a computational modeling approach called information flow (IF) analysis62,63. To perform IF analysis from the EGFR at the cell membrane to nuclear transcription factors (TFs), the two EGFRNets were first supplemented with publicly available prey–prey interactions from InnateDB.com64 and 122 additional nodes that are known to be involved in EGFR signaling18 but were not chosen as bait proteins in our AP-MS experiments (Supplementary Data 14). These networks are referred to as the HCT116IFANET and HKE3IFANET. Information flow analysis was implemented using the CytoITMprobe software (damping factor = 0.85; channel model selected)41, selecting EGFR as the source node of signaling and 19 downstream TFs (Supplementary Data 14) as the sinks for the information flow. Information flow scores were determined by measuring how much information flows through each node in the HCT116IFANET and HKE3IFANET networks. See Supplementary Methods, section 26, for a detailed description.

Gene ontology, pathway, and transcription factor binding site analyses

Gene ontology (GO) and pathway analyses were performed using InnateDB.com64 regarding GO terms or pathways with an FDR < 0.05 as significantly enriched. Transcription factor binding site analysis was undertaken using the findMotifs.pl program in HOMER v4.865, with the human hg38 promoter set in order to identify enriched motifs. See Supplementary Methods, section 27, for a detailed description.

Analysis of CRC patient data

Survival data of 629 CRC patients were obtained from the TCGA50, and correlated with alterations in genes encoding either the top 20 most rewired or the bottom 20 least rewired bait proteins. The alterations included were mutations, copy number changes, mRNA expression changes, and protein expression changes. As an additional control we selected a set of 36 baits that accounted for the same number of interactions as the top 20 baits in the network. The Kaplan–Meier curves were plotted using PRISM 7.0.3. See Supplementary Methods, section 27, for a detailed description.

Western blotting

Cells were lysed in 1% NP40, 20 mM Tris-HCl pH 7.5, 150 mM NaCl, 1 mM MgCl2) supplemented with protease inhibitor cocktail (Roche) and phosphatase inhibitors (2 mM sodium orthovanadate, 10 mM sodium fluoride and 10 mM β-glycerophosphate; all from Sigma-Aldrich) for 10 min at 4 °C. Lysates were cleared by centrifugation at 20,000 × g for 10 min, and adjusted to equal protein concentrations. Proteins were separated by sodium-dodecylsulfate polyacrylamide gel electrophoresis (SDS-PAGE) and transferred to polyvinylidene difluoride (PVDF) membranes. Blots were incubated with the respective antibodies and developed using Enhanced Chemiluminescence (ECL; Thermo Fisher) according to the manufacturer’s instructions. Blots were quantified using the Image J software and phospho-specific antibody signals were normalized to the total abundance of the respective proteins.

Luciferase assays

The transcription factor response element activity was assessed with luciferase constructs bearing response elements for STAT1 (4xGAS response element; Stratagene, #219091–51), STAT1/STAT2 (IRSE/interferon alpha response element66) and STAT3 (4xm67 response element67). A CMV-β-gal plasmid was co-transfected as a control of transfection efficiency. Forty-eight hours post transfection cells were stimulated with 10 nm human EGF (Roche; #11376454001) for 5 h before luciferase and β-gal activity were measured using luciferase assay (Promega, #E4030) and β-galactosidase assay kits (Promega, #E2000). Luciferase activity was normalized againstβ-gal activity to correct for transfection efficiency.

Transcriptional profiling

HCT116 (mtKRASHi) and HKE3 (mtKRASLo) cells were serum starved for 18 h before stimulation with TGF-α (0.01 µg/ml, Abcam) for 0, 15, 30, 60, 90, and 120 min. RNA was extracted using the TRIzol reagent (Thermo Fisher Scientific) from three biological replicates at each time point. RNAseq was performed with an Illumina HiSeq 2500 machine using a v2 High Output 100 cycle Kit (1 × 100 bp SR). Data were processed as described in the Supplementary Methods.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.