Magnitude of Ubiquitination Determines the Fate of Epidermal Growth Factor Receptor Upon Ligand Stimulation

Receptor tyrosine kinases (RTK) bind growth factors and are critical for cell proliferation and differentiation. Their dysregulation leads to a loss of growth control, often resulting in cancer. Epidermal growth factor receptor (EGFR) is the prototypic RTK and can bind several ligands exhibiting distinct mitogenic potentials. Whereas the phosphorylation on individual EGFR sites and their roles for downstream signaling have been extensively studied, less is known about ligand-specific ubiquitination events on EGFR, which are crucial for signal attenuation and termination. We used a proteomicsbased workflow for absolute quantitation combined with mathematical modeling to unveil potentially decisive ubiquitination events on EGFR from the first 30 seconds to 15 minutes of stimulation. Four ligands were used for stimulation: epidermal growth factor (EGF), heparin-binding-EGF like growth factor, transforming growth factor-a and epiregulin. Whereas only little differences in the order of individual ubiquitination sites were observed, the overall amount of modified receptor differed depending on the used ligand, indicating that absolute magnitude of EGFR ubiquitination, and not distinctly regulated ubiquitination sites, is a major determinant for signal attenuation and the subsequent cellular outcomes. 2021 The Author(s). Published by Elsevier Ltd. This is an open access article under the CCBY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Epidermal growth factor receptor (EGFR) is a member of the ErbB family of receptor tyrosine kinases (RTKs) and is involved in the regulation of vital cellular processes most notably cell proliferation and differentiation. 1,2 Its dysregulation, e.g. as a result of overexpression or activating mutations, is involved in the genesis of many malignancies. 3,4 So far, seven different ligands of EGFR have been identified: EGF, heparin-binding EGFlike growth factor (HB-EGF), epiregulin (EPI), transforming growth factor-a (TGFa), betacellulin, amphiregulin, and epigen. 5 Binding of a ligand to the extracellular parts of EGFR triggers its homo-or heterodimerization with other ErbB receptors and subsequent intracellular autophosphorylation, 6 which results in increased enzymatic activity and initiation of a signaling cascade at the plasma membrane. 1 Simultaneous with signal initiation, its attenuation and termination are triggered as well. The endocytic removal of active receptor-ligand complexes from the cell surface followed by their lysosomal degradation plays a major role in this process. Receptor internalization has been shown to happen via clathrin-mediated and clathrinindependent endocytosis pathways and is influenced by ligand concentrations and the employed experimental system (for review see 7 ). A common molecular process involved in the majority of proposed pathways is the covalent conjugation of ubiquitin (monoubiquitination) or a chain of ubiquitins (polyubiquitination) to the receptors. [8][9][10][11][12] Previous studies have demonstrated both K63type polyubiquitination 13 as well as EGFR monoubiquitination on multiple sites, 14,15 to be responsible for endocytosis and degradation of the receptor.
A phenomenon, which is only understood in few cases, is the observation that different ligands of the same RTK can cause distinct cellular responses encoded by ligand-specific signaling strength and duration. 16 For EGFR it has been shown that EGF mainly induces receptor degradation, whereas TGFa leads to pronounced receptor recycling, 17,18 which probably contributes to the observed stronger mitogenic activity of TGFa. 19 Using an integrated multilayered proteomics approach, the differential effects of EGF and TGFa were linked to the differential phosphorylation of a down-stream signal transducer, the vesicle trafficking protein RAB7, and differential recruitment of the EGFR interacting protein RAB11FIP1. However, an enigma remains the primary, initial molecular event leading to the differential downstream effects, which are observable after several minutes of ligand stimulation. At these time points more than 50% of cell surface EGFR are already internalized. 19 The most proximal event after ligand binding is receptor dimerization and the occurring of posttranslational modifications, and it can be assumed that binding of distinct ligands leads already to differential effects on the receptor itself. Two scenarios are conceivable: due to steric differences, different ligands may lead either to (i) modifications of distinct sites in the cytosolic domains of the receptor pairs, or (ii) distinct absolute abundances of modifications on common sites, which are encoding the specific downstream events. To determine if specific ubiquitination sites or the extent of ubiquitination is critical for distinct receptor internalization routes followed by termination of signaling, we absolutely quantified the ubiquitination events within the first 15 min of EGFR stimulation using four ligands with different mitogenic potentials: EGF, HB-EGF, TGFa and EPI.

Results
The growth factors EGF, HB-EGF, TGFa and EPI exhibit differential mitogenic potential As model system, we chose the human cell line Hep2, which has been used in the past to describe distinct, ligand-dependent internalization routes of EGFR. 17 As EGFR internalization is growth factor (GF) concentration dependent 9 and to prevent potential differences due to partial receptor activation, we decided to choose GF concentrations which induce a robust, comparable receptor activation/phosphorylation, although these concentrations are higher than the levels observed in vivo. We tested three concentrations per GF, performed EGFR immunoprecipitations following 6 min of stimulation, and analyzed the levels of receptor tyrosine autophosphorylation by western blot, representing the most proximal readout for receptor activation (Supplemental Figure S1(a)). 20 In the case of EGF, HB-EGF and TGFa, 150 ng/ml resulted in a maximal, robust modification of the receptor, whereas 500 ng/ml of EPI were needed to achieve a similar extent of EGFR phosphorylation. Such GF concentration dependencies were described previously for these cells. 17 Next to western blot analyses of global receptor phosphotyrosine levels, we also used synthetic peptides to absolutely quantify specific autophosphorylation events by mass spectrometry (MS). In agreement with previously published data, 19 we detected the maximum autophosphorylation, corresponding to ca. 60% of receptor molecules, after 2-6 min of stimulation ( Figure 1(a), Supplemental Table S1). Although we used more than three times the concentration of the other three GFs, EPI appeared to be less potent leading to phosphorylation events on only 40% of EGFR molecules. Otherwise, the kinetics of phosphorylation events were comparable and exhibited only minor GF-dependent differences.
To study the phenotypic consequences of GF treatments, we stimulated Hep2 cells with respective GF concentrations and analyzed their effects on cell proliferation and receptor recycling. Corroborating published results, TGFa and EPI were the strongest mitogens, followed by EGF and HB-EGF (Figure 1(b), Supplemental Figure S1(b)). [17][18][19] Furthermore, the mitogenic potentials of the GFs correlated well with receptor recycling/surface presentation: whereas 48% and 67% of total EGFR were present at the cell surface 45 min after stimulation with TGFa and EPI, respectively, this was the case for merely 23% and 11% of EGFR molecules after EGF and HB-EGF treatment, respectively (Figure 1(c), Supplemental Table S2).
Thus, although the different GF lead to similar extents and kinetics of receptor autophosphorylation, they procure different phenotypic outcomes. As receptor recycling and trafficking are linked to its ubiquitination, 10,14,21 we decided to study GF-induced receptor ubiquitination dynamics in detail to define potential GF-specific effects that could lead to the observed phenotypes.

Growth factors induce differential extent of EGFR ubiquitination
In MS-based proteomics experiments the digestion of proteins by trypsin, which generates peptides readily analyzable by liquid chromatography (LC)-MS/MS, is a standard working procedure. Tryptic digestion of ubiquitinated proteins yields branched modified peptides with the C-terminal di-Glycine (-GG) remnant of ubiquitin being attached to the lysine residue carrying the modification, which is used in MS analyses to identify substrate ubiquitination sites. In this manner, 6 ubiquitination sites on human EGFR have been initially deposited into UniProt, 21 all of them situated on the surface of the kinase domain (Figure 2(a)). Our aim was to absolutely quantify stimulus-specific ubiquitination kinetics on EGFR within the first 15 min of receptor stimulation. To achieve this goal we used synthetic, -GG-modified peptides containing the known ubiquitination sites K716, K737, K754, K867, K929, K970. 13,22 In addition, we chose 6 peptides containing potential ubiquitination sites, which were identified as being ubiquitinated in different MS screens [23][24][25][26] : K713, K757, K846 lying within the kinase domain, K708 just in front and K1061 and K1188 just after the kinase domain ( Figure 2(a), Table 1). We then utilized a "reverse AQUA" method 20 in combination with SILAC labeling 22 to accurately determine the kinetics of these ubiquitination sites on EGFR following stimulation with the four GFs ( Figure 2(b)). The original AQUA approach relies on a heavy, stable isotopelabeled, synthetic peptide for absolute quantitation of an endogenous peptide of interest. The synthetic AQUA and the native peptide have identical sequences but differ in mass, which makes them easily distinguishable in mass spectrometers. Spiking known amount of the AQUA peptide into a complex biological sample then allows to precisely measure the quantity of the native peptide in that sample by MS. 20 In the "reverse AQUA" method here, we used regular synthetic peptides (not labeled with isotopes), and we spiked those in samples originating from 'medium-heavy' and 'heavy' SILAC-labeled cells. 22 In this way, we could do absolute quantitation of our peptides of interest in two biological samples simultaneously. Hep2 cells were SILAC labeled by medium-heavy and heavy arginine and lysine variants and five SILAC experiments for each GF were performed, yielding altogether biological replicas for 0 s, 30 s, 2 min, 6 min, and 15 min time-points of ligand stimulation ( Figure 2(b)). Lysates of each experiment were combined and EGFR purified by immuno-affinity purification using anti-EGFR antibodies followed by SDS-PAGE and in-gel digestion of the upper gel region containing EGFR and the posttranslationally modified variants of the receptor (Supplemental Figure S2). Synthetic, unlabeled peptides were spiked in known amounts to the tryptic peptide mixtures allowing the combination of data and construction of five time-point kinetics. In the chosen experimental approach, absolute quantitation of ubiquitination events is based on comparing the intensities of the endogenous and the corresponding synthetic peptides in the same spectra. However, if a fraction of the endogenous ubiquitinated peptide carries an additional modification, e.g. phosphorylation, such doublymodified peptide will have a different mass and chromatographic property and will compromise the accuracy of quantitation. Thus, recorded ubiquitination dynamics might be obscured by additional phosphorylation events on the respective endogenous peptides. This could be the case for example with the tryptic peptide bearing K867, which also contains the known phosphorylated amino acid residue Tyr869. To rule out this potential error we stimulated SILAClabeled cells with EGF, treated one sample with lambda phosphatase, and quantified ubiquitination events before and after treatment (Supplemental Figure S3). The removing of phosphate groups had no detectable influence on the intensities of ubiquitinated peptides.
Using 150 ng/ml of TGFa, EGF and HB-EGF and 500 ng/ml EPI we quantified the ubiquitination events on EGFR up to 15 min of stimulation ( Figure 3, Supplemental Table S1). As expected, compared to the tyrosine phosphorylation sites (Figure 1(a)), ubiquitination sites responded substantially slower, with high levels of ubiquitination being observable at the 6 and 15 min time points. Also, the extent of ubiquitination was lower compared to phosphorylation, the peak amount of ubiquitination being detected on K737 after HB-EGF treatment for 6 min (25%, Figure 3). Comparing different GFs and time points, ubiquitination on K708, K846, and K1188 appeared to be unresponsive in all cases. The relative order of modification by ubiquitination remained the same across the responsive sites regardless of the used GF: K737 showing the highest ubiquitination, followed by K716, K754 and K867. The most pronounced difference between GF treatments were the extents of ubiquitination. HB-EGF, having the lowest mitogenic and recycling activity, exhibited the strongest ubiquitination, whereas TGFa and EPI, having the highest mitogenic activity, exhibited the weakest ubiquitination levels on all sites.

All growth factors induce K63 ubiquitination
Next to ubiquitination sites on EGFR itself, we were also interested in the K63 ubiquitin chain linkage, as this type of polyubiquitination was shown previously to occur on EGFR upon ligand stimulation. 13,22 For that purpose, in addition to all spiked AQUA peptides corresponding to the different ubiquitination sites on EGFR, we added a -GG-modified synthetic peptide matching the characteristic peptide derived from the tryptic digest of K63 polyubiquitin chains (Supplemental Table S1). Furthermore, we also monitored the K6, K11 and K48 ubiquitin chains, all by the aid of synthetic peptides mimicking the unique tryptic product for each of the corresponding chain types 27 (Supplemental Table S1). All synthetic peptides were spiked Table 1 Synthetic peptides used for the absolute quantitation of EGFR ubiquitination sites. Synthetic peptides were spiked into eluates of in-gel digests for mass spectrometric analyses. Ubiquitinated lysine residues are marked with an asterisk. Carbamidomethylated cysteines are indicated as C # . together in the samples derived from the 'mediumheavy' and 'heavy' SILAC-labeled cells stimulated with the GFs (Figure 2(b)). In the EGFR immunoaffinity purifications, the fractions of K6 and K11 chains remained unaffected by the GF treatments, K48 chains showed small decreases, whereas the K63-linked peptides increased over time with each of the ligands (Figure 4(a)). Notably, HB-EGF and EGF stimulations triggered stronger K63 chain formations than TGFa and EPI, accounting for about 50% of the total EGFR ubiquitination at the 6 min time point (Figure 4(a)). Similarly, quantifying the total level of ubiquitin detected in the EGFR immuno-affinity purifications revealed that HB-EGF led to the most pronounced increase, followed by EGF, TGFa and EPI, with the latter two inducing less than half of the total ubiquitination detected after 6 min of HB-EGF stimulation (Figure 4(b), Supplemental Table S1). Thus, it appeared that not the ubiquitination of specific sites, but rather the magnitude of ubiquitination of common sites determined the routes of EGFR trafficking and by this the cellular phenotype. To assess whether ubiquitination magnitude could indeed code for the observed phenotypes, we turned to mathematical modeling.
Mathematical modeling highlights that quantitative differences in common signals may lead to discrete trafficking routes A mathematical model was built by a system of ordinary differential equations (ODEs) describing the relevant receptor reactions based on massaction kinetics (Supplemental Data). It was developed to describe the main effects initiated by ligand stimulation of the cells, testing whether a difference in ligand affinity was sufficient to describe phenotypic differences. The model was kept as simple as possible to avoid over-fitting. The receptor states were classified using the receptor properties of ligand-association, phosphorylation, ubiquitination and localization status, resulting in a total of 2 4 = 16 states (Supplemental Figure S4). The production of new receptors was neglected and the model should be interpreted as the deviation from the equilibrium state of the cells. The model was calibrated using the measurements of the phosphorylated, ubiquitinated and membrane bound receptor fractions, denoted as EGFR p , EGFR u and EGFR m . The predicted trajectories and data points are shown in Figure 5(a) (solid line). In general, the model was able to capture the main differences of the data. Small deviations, however, are visible, especially for the peak of the ubiquitination for the HB-EGF data. This peak could also not be explained by a HB-EGF specific behavior by allowing the rate constants to differ for each ligand (see dashed line Figure 5(a)). Furthermore, the following modifications of the model were tested. The environment of the internalized receptors differs from the environment of membrane bound receptors. This could lead to a change in the rate constants for ligand dissociation, (de)phosphorylation and (de) ubiquitination, implemented as first modification. Second, the distinction of recycling versus degradation of internal receptors based purely on the ubiquitination on-or off-state of the internal receptors might be a too strong assumption. Therefore, two more rates constants were introduced, k r2 for the recycling of ubiquitinated receptors and k d2 for the degradation of nonubiquitinated receptors (see dotted lines in Figure 5(a)). These modifications allowed for a slightly better description of the TGFa data but did not lead to a better fit of the HB-EGF data.
The correct description of receptor trafficking processes with the EGFR m data as readout can be interpreted as a phenotypic validation of the model ( Figure 5(a), lowest panels). In addition, we also used the cell proliferation measurements and tested if the model, calibrated on receptor data for up to 45 minutes of stimulation, could describe the ligand dependent proliferation of the cells (Supplemental Data). The model could recapitulate the differences observed in liganddependent proliferation data as well ( Figure 5(b)), further indicating that differences in signal strength caused by a different ligand affinity correlates well with the observed ligand-specific cellular phenotypes.

Discussion
Due to its high biological relevance and involvement in numerous human malignancies, EGFR signaling is one of the best studied signal transduction pathways, also by quantitative MSbased proteomic studies. [28][29][30] While a lot is known about the role of individual phosphorylation sites on EGFR, site-specific ubiquitination events involved in downregulation of signaling are far less understood. Using synthetic peptides, we quantified the absolute levels of individual ubiquitination events on the receptor and detected only minor qualitative ligand-specific differences arguing that ubiquitination strength, rather than site-specificity, is primarily responsible for the differential downstream effects on EGFR trafficking and degradation.
In addition to EGFR homodimerization, some ligands can also trigger EGFR heterodimerization with other members of the ERBB receptor family, which could potentially contribute to the downstream phenotypic effects. However, this is not likely to be the case for the Hep2 cells as we only detected negligible traces of ERBB2 and no measurable amounts of the ERBB3 and ERBB4 in the EGFR immunoprecipitations from the cells stimulated with these growth factors. In other cell types or cell line models heterodimerization may nevertheless play an important role.
As absolute quantitation by MS depends on detecting the exact endogenous counterpart of the spiked synthetic peptide, respective endogenous peptides carrying additional PTMs are not included in the calculations and will obscure quantitation results. In our case, endogenous ubiquitinated peptides additionally carrying phosphate groups would lead to an underestimation of the amount of a specific ubiquitinated site. As the used ubiquitinated peptides include potential phosphorylation sites, 31 we were concerned that these potentially doubly modified peptides would interfere with our experiments. To address this question, we treated one sample with lambda phosphatase removing all attached phosphate groups and compared the intensities of ubiquitinated peptides to untreated control samples. In all instances no difference before and after lambda phosphatase treatment was observed, supporting our experimental approach. In addition, these results indicate that neighboring phospho-and ubiquitin-sites might not be modified simultaneously, most likely due to sterical interference, either of the PTM itself, the modifying enzymes or of interacting proteins binding to the respective sites. In this regard, EGFR has also been reported to bear O-N-Acetyl glucosamine (O-GlcNAc) modifications that could have affected the accuracy of quantitation as well, if occurring at positions nearby the examined sites of ubiquitination. 32 However, it appeared that it was not the case here as only two O-GlcNAc sites on the EGFR were identified when considering this PTM as a possible variable modification as well, and both positions (T354 and S380) were located on the extracellular part of the receptor.
EGFR has multiple lysine residues, which may all get ubiquitinated in specific experimental settings. A variant missing 21 lysine residues in its kinase domain alone was used to study EGFR internalization and was found to still exhibit liganddependent, residual internalization rates, 33 which could finally be linked to cryptic/less abundant ubiquitination sites in the receptor itself. 9 Thus, by following 12 ubiquitination sites we do not make the claim to cover comprehensively all possibilities. However, we have very likely covered the major ligand-responsive sites as we did not identify other ubiquitination sites on EGFR by shot-gun proteomics in these samples. Any ligand-induced ubiquitinations, which we might have missed, probably affect only a minor fraction of receptor molecules and will not be able to account for the observed phenotypic differences.
We were surprised to find the same relative order of ubiquitination on individual sites comparing four different GF. The overall dynamics of observed ubiquitination events were also quite similar, whereas the amount of modified receptor varied. Thus, in contrast to phosphorylation events, which have been shown to vary qualitatively between different ligands in other RTKs, 16 ubiquitination events in EGFR signaling appear to transfer information by the magnitude of modification. This would also correlate with the observation that nonfunctional individual or combinatorial mutations of EGFR ubiquitination sites have less pronounced biological effects. 34 Finally, using the generated quantitative time-resolved data of EGFR phosphorylation, ubiquitination and ubiquitin chainlinkages, we developed a mathematical model based on a system of ODEs, which support analyses of the dynamic behavior of cellular systems. ODE-based models have a high predictive power allowing non-intuitive insights into complex biological systems. 35 In our case, we used the model in a reductionist-fashion and asked if the magnitude of ubiquitination in contrast to ubiquitination events targeting different sites within receptor molecules could be a determinant for ligand-specific receptor degradation rates resulting in altered cell proliferation. Indeed, the tested model variants argue for the hypothesis that ubiquitination strength can be a determinant for the observed GF-dependent phenotypic differences. Nevertheless, the model also indicated that additional research is necessary as we could not recapitulate in detail the HB-EGF induced ubiquitination dynamics. Although it is not clear what are the factors controlling the level of receptor ubiquitination, a key feature may be the strength of EGFR dimerization triggered by the different ligands. It has already been shown that EPI stimulation induces weaker and shorter-lived EGFR dimers compared to EGF and, in general, EGFR dimer stability can vary significantly depending on the ligand it is bound to. 36 According to our model, however, different strengths of receptor interactions by themselves, i.e. different on and off rates of the ligands, hold not enough information to explain all observations.
In conclusion, we have generated absolute quantitative data of ligand-specific EGFR ubiquitination events and in combination with mathematical modeling, we could show that the overall amount of modified receptor (the magnitude of receptor ubiquitination) appears to be a major determinant for the different endocytosis phenotypes and subsequent cellular outcomes induced by the growth factors.

Materials and methods
Experimental design and statistical rationale Experiments for assessing cell proliferation were performed in biological triplicates for each condition and the entire experiment was repeated once more. The mathematical modeling of the data is described in details in a separate section.

Lysis of cells and immunoprecipitation of EGFR
After stimulation cells were thawed and lysed in ice-cold lysis buffer containing 150 mM NaCl, 1 mM EDTA, 50 mM Tris, pH 8.0, 0.5% NP-40, 0.5% Triton X-100, 0.25% sodium deoxycholate, 1 mM sodium orthovanadate, 5 mM NaF, 5 mM bglycerophosphate, and protease inhibitors (Complete tablets, Roche Diagnostics). Prior to cell lysis chloroacetamide was added to the lysis buffer (30 mM) to inhibit deubiquitination enzymes. For the MS-based experiments lysates with different labeled cells were combined as described in the preceding paragraph and centrifuged at 14'000 g. Supernatants were used for immunoprecipitation with 4 lg of pre-bound Protein G Sepharose 4 Fast Flow (GE Healthcare) anti-EFGR antibody per 15 cm dish (mAb Ab-11 clone 199.12 (Invitrogen), mAb 528 sc-120 (Santa Cruz Biotechnology)) for 5 h at 4°C followed by washing the beads five times with lysis buffer.

MS sample preparation and analysis
Beads were incubated in SDS sample buffer for 10 min at 90°C, proteins eluted and resolved on Novex 4-12% Bis-TRIS gradient gels using the MES buffer system (Invitrogen) followed by staining of the gel (Colloidal Blue Staining Kit, Invitrogen). The region of the gel containing unmodified and modified EGFR was excised and proteins were subjected to in-gel digestion essentially as described 37 with modifications of using chloroacetamide for alkylation instead of iodoacetamide 38 and a mix of LysC enzyme (12.5 ng/ml) and Trypsin (12.5 ng/ml). The mastermix of synthetic peptides (containing 1 pmol of each synthetic peptide) was spiked into samples prior to the extraction step of the in-gel digestion procedure. The extracted peptide mixtures were concentrated and desalted using STAGE tips. Prior to MS analysis samples were dissolved in a solvent of 0.1% TFA containing 0.01% H 2 O 2 as it was used here 39 and injected into a 24 cm fused silica column with an inner diameter of 75 mm packed in-house with C18 resin (3 mm beads, Reprosil, Dr. Maisch GmbH) for reverse-phase chromatography using an EASY-nLC system (Thermo Fisher Scientific) that was connected on-line to a Q Exactive mass spectrometer (Thermo Fisher Scientific) equipped with a nano-electrospray ion source (Thermo Fisher Scientific). Peptides were loaded in solvent A (0.5 % acetic acid) and eluted by applying a 120 min gradi-ent of solvent B (80% ACN, 0.5% acetic acid). The Q Exactive mass spectrometer was operated in positive polarity mode with a capillary temperature of 275°C. Full MS survey scan resolution was set to 70'000 with automatic gain control (AGC) target value of 1e6 for a scan range of 300-1750 m/z and maximum ion injection time (IT) of 120 ms. A data-dependent method was used for acquisition: the top 12 most intense ions were fragmented by higher-energy collisional dissociation (HCD) with a normalized collisional energy (NCE) of 25 eV. Precursor ions with charge states 1, 8 and higher were excluded from selection. MS/MS scans were performed with a resolution of 35'000, maximum IT of 124 ms and an ion target value of 5e5, scan range of 200 to 2000 m/z, 1.2 m/z isolation window. Repeat sequencing of peptides was prevented by setting the dynamic exclusion window to 45 seconds.
A single peak-list was generated from raw files using DTA supercharge v.

Lambda phosphatase treatment
Hep2 cells were SILAC labeled as described above and stimulated with EGF for 6 minutes. Lysates from differently labeled cells were handled separately for further treatment. EGFR was immunoprecipitated as described above, sepharose beads (50 ml) were washed with lysis buffer without EDTA, phosphatase and protease inhibitors. The IP of heavy labelled cells (Arg10, Lys8) was incubated with 1000 units (AU) of recombinant Lambda Protein Phosphatase (k-PPase, Millipore) in the presence of 100 ml of k-PPase buffer for 40 minutes at 30°C, whereas the IP of medium labelled cells (Arg6, Lys4) was mock-treated for the same time at 30°C. After adding of SDS sample buffer k-PPaseand mocktreated samples were combined and resolved on SDS-gel for further analysis as described above.

Proliferation assay
The analysis of cell proliferation in response to four different growth factors was performed as described 41 . Briefly, we used cells counting and seeded 33,000 Hep2 cells per well for 8 hours. Cells were then washed 2 times in serum-free medium and stimulated with GFs in serum-free medium for 16 and 24 hours. After incubation cells were detached from plates and counted using NucleoCounter Ò NC-100 TM (Chemometec) in accord to manufacturer's instructions. We performed three independent biological replicas for each time-point and calculated standard deviation for each time-point measurement. Experiments were performed twice.
Measuring EGFR on the cell surface Hep2 were grown in 6 well plates to 90% confluence and starved overnight. Cells were incubated with GF for different time points (0 (unstimulated control), 15 and 45 minutes), washed with 3 ml of ice-cold PBS twice and plates were kept on ice. The ice-cold PBS solution of antibodies (1 mg/ml of A11 (Invitrogen) in combination with 1 mg/ml of 528 from Santa Cruz), was applied (1 ml to 1 well) and incubated for 1.5 h at 4°C gently shaking. Cells were washed with 3 ml of ice-cold PBS twice and lysed with 700 ml of the lysis buffer as described above. The lysates were cleared by centrifugation, supernatants were added to a 10 ml mix of Protein A-G Sepharose beads. The IPs were incubated for 2 h at 4°C and washed 5 times with ice-cold PBS. Proteins were eluted from beads with 50 ml of 8 M Guanidine HCl incubating the mix at 95°C for 5 min. After elution proteins were subjected to in-solution digest with the LysC and Trypsin overnight as described in, 42 purified with the STAGE-tips and run on a LTQ Orbitrap Velos Pro MS-instrument (Thermo Scientific) using Parallel Reaction Monitoring (PRM).

PRM analysis
All PRM analyses were performed similar to 43 with minor modifications. An LTQ Orbitrap Velos Pro MS was coupled with the same LC system described before for Data Dependent Analysis (DDA). Precursor m/z of 9 proteotypic target peptides for EGFR, (Supplemental Table S2) have been selected from previous DDA runs and added to the instrument inclusion list. Targeted MS/MS spectra were acquired in the linear ion trap using a global unscheduled inclusion lists where every peptide was fragmented in CID mode (CE = 35). Isolation window and activation time were set respectively to 2 Da and 10 ms. An acetonitrile gradient (from 3 to 45% in 60 min) was applied to obtain an optimal number of data point (!10) along the peak elution profile for the quantification. PRM data analysis was carried out on Skyline TM 3.6.0 software. 44 Spectral libraries were built using Max-Quant ms/ms search files. Resolving power and mass analyzer were set respectively to 0.7 m/Z and quadrupole ion trap (QIT). Area under the curve (AUC) values relative to the 5 most intense product ions for each peptide were exported and plotted in a matrix, which was used for the quantitative comparison (Supplemental Table S2).

Mathematical modeling
The mathematical model is described by a set of reactions that are mathematically formulated as ordinary differential equations (ODEs), following mass-action kinetics. The full set of ODEs is given in the Supplemental Data. The states of the model are linked to the data by an observation function. The observables EGFR p , EGFR u and EGFR m are defined as the sum of all phosphorylated, ubiquitinated or membrane bound receptor states, devided by the sum off all receptor states in case of EGFR p and EGFR u . In case of the phosphorylation and ubiquitination data, proportionality factors were implemented to scale the model states to the data: EGFR data p = sp•EGFR p and EGFR data u = su•EGFR u . A model prediction is computed for a set of parameters by solving the ODEs numerically. The parameters contain values for the rate constants associated with the reactions, initial values for the receptor states and ligand and the proportionality factors.
The parameters are determined from the data by the maximum-likelihood approach, using the where y(t i , h) is the prediction of observable i at time point t i given the parameters h and y D i are the data points with uncertainty r i at time point t i . Parameter estimates are obtained by minimizing Eq. (1) with respect to h. The profile-likelihood method as presented in 46 is used to evaluate the identifiability of the parameters.
The model-based analysis, including the parameter estimation and identifiability analysis, has been performed using the dMod 47 package for R. The package is available on the Comprehensive R Archive Network (CRAN).

Data Availability Statement
All mass spectrometric data generated in this study have been submitted to the ProteomeXchange Consortium via the PRIDE partner repository with the data set identifier PXD019621.