Time-resolved Phosphoproteome Analysis of Paradoxical RAF Activation Reveals Novel Targets of ERK *

Small molecules targeting aberrant RAF activity, like vemurafenib (PLX4032), are highly effective against cancers harboring the V600E BRAF mutation and are now approved for clinical use against metastatic melanoma. However, in tissues showing elevated RAS activity and in RAS mutant tumors, these inhibitors stimulate RAF dimerization, resulting in inhibitor resistance and downstream “paradoxical” ERK activation. To understand the global signaling response of cancer cells to RAF inhibitors, we profiled the temporal changes of the phosphoproteome of two colon cancer cell lines (Colo205 and HCT116) that respond differently to vemurafenib. Comprehensive data mining and filtering identified a total of 37,910 phosphorylation sites, 660 of which were dynamically modulated upon treatment with vemurafenib. We established that 83% of these dynamic phosphorylation sites were modulated in accordance with the phospho-ERK profile of the two cell lines. Accordingly, kinase substrate prediction algorithms linked most of these dynamic sites to direct ERK1/2-mediated phosphorylation, supporting a low off-target rate for vemurafenib. Functional classification of target proteins indicated the enrichment of known (nuclear pore, transcription factors, and RAS-RTK signaling) and novel (Rho GTPases signaling and actin cytoskeleton) ERK-controlled functions. Our phosphoproteomic data combined with experimental validation established novel dynamic connections between ERK signaling and the transcriptional regulators TEAD3 (Hippo pathway), MKL1, and MKL2 (Rho serum-response elements pathway). We also confirm that an ERK-docking site found in MKL1 is directly antagonized by overlapping actin binding, defining a novel mechanism of actin-modulated phosphorylation. Altogether, time-resolved phosphoproteomics further documented vemurafenib selectivity and identified novel ERK downstream substrates.

Increased RAS-ERK signaling is critically connected with tumor initiation and progression (1). Canonical pathway activation is initiated at the plasma membrane where ligand bound to receptor tyrosine kinases (RTKs) 1 transduce the signal to enzymes that activate RAS (NRAS, KRAS, or HRAS) by favoring its GTP loading (2). Downstream of RAS, the RAF family of Ser/Thr kinases comprises three different homologs, ARAF, BRAF, and CRAF, involved in the regulation of the RAS-ERK signaling cascade (1,3). In quiescent cells, RAF is kept inactive in the cytoplasm in an auto-inhibited complex with 14-3-3 proteins. GTP-bound RAS directly associates with the RAF RAS-binding domain, which recruits RAF to the plasma membrane and relieves its auto-inhibition (3,4). Importantly, RAF activation also relies on the formation of homo-and heterodimers that allosterically control kinase activity (5). Upon activation, RAF signals downstream to MEK1/2 (MAP2K1/2), which in turn phosphorylate ERK1/2 (MAPK1/3) proteins leading to the phosphorylation of hundreds of cellular substrates on consensus PX(S/T)P motifs (6).
RAS and RAF oncogenic mutations render them constitutively active resulting in aberrant signaling to downstream effectors. This induces a sustained activation of downstream MEK and ERK in a growth factor-independent fashion, a hallmark of cancer cell signaling. The most frequent genetic lesions in RAS genes, the mutation of Gly-12, Gly-13, and Gln-61 residues, lock RAS proteins in the GTP-bound active state (2). Similarly, BRAF mutations like V600E occur in the kinase catalytic domain and result in a conformational change rendering BRAF constitutively active (7,8). RAS mutations are one of the leading causes of cancer irrespective of tissue type (9), whereas BRAF V600E mutations are also found across several tumor types but are especially frequent in melanoma (70%) and thyroid cancers (45%) (10 -12).
In view of the frequent association of RAF with cancer signaling, several drug discovery programs have aimed at blocking its kinase activity (13). Pharmacological inhibitors targeting oncogenic RAF such as the Food and Drug Administration-approved PLX4032 (vemurafenib, Vectibixா, or Erbituxா) can abrogate the activity of the V600E mutant BRAF (14 -16). However, wild type (WT) RAF tumors bearing RAS mutations or showing elevated RTK signaling give rise to aberrant activation of RAF kinases upon RAF inhibitor treatment (17)(18)(19). This undesired effect is caused by a RAS-dependent increase in RAF dimerization upon inhibitor binding, which leads to the so-called "paradoxical" activation of ERK signaling. The molecular basis of this phenomenon is the allosteric activation of drug-free RAF by drug-bound protomers across induced dimers (17)(18)(19)(20). Several cancer types therefore display poor response to RAF-specific inhibitors. This limits the use of ATP-competitive inhibitors like vemurafenib (PLX4032) to the treatment of melanomas carrying the BRAF V600E mutation that have reduced RAS signaling (21). Cancer cell lines derived from RAS and RAF mutant tumors thus behave in an opposite manner when exposed to RAF inhibitors. Although the paradoxical effect of RAF inhibitors is undesired in clinical settings, it represents a powerful tool to selectively induce pathway activation without the pleiotropic effects typically associated with RTK or phorbol ester stimulations.
We recently showed that dynamic profiling of phosphoproteomes enabled the identification of novel downstream effectors of well studied signaling pathways with high confidence (22,23). This type of proteomics analysis facilitated the characterization of biological responses to a variety of stimuli like osmotic shock and thermal stress (22)(23)(24). Phosphoproteome dynamics provides valuable information on the shape and extent of a biological response to unravel novel and unexpected kinase-substrate interactions (22). In the specific context of the RAS-ERK signaling cascade, the identity of ERK substrates and the timing of pathway activation determine what biological response will be elicited (proliferation, differentiation, and survival) (25). Systems approaches that include treatment kinetics are therefore needed to derive an exhaustive repertoire of ERK substrates.
Although phosphoproteome analysis of RAS-ERK pathway blockade has recently been published (26), the proteomewide consequences of the paradoxical pathway activation have not yet been investigated. We posit that a dynamic understanding of this phenomenon could help decipher the complex response to RAF inhibitors and yield new insights on RAS-ERK pathway signaling. In this study, we report the application of high resolution temporal profiling coupled with metabolic labeling (27) to quantify changes in the phospho-proteome of HCT116 and Colo205 colon cancer cells following vemurafenib treatment. Altogether, this approach enabled the identification of 37,910 phosphorylation sites, of which 660 were dynamically regulated. We used this rich dataset to understand the complex response to RAF inhibitors and to identify additional signaling nodes that act downstream of ERK.

EXPERIMENTAL PROCEDURES
Alphascreenா SureFireா p-ERK1/2 Assay-The pERK Alphascreenா SureFireா (PerkinElmer Life Sciences) was used following the manufacturer's instructions for the standard two-plate assay protocol. For the time course experiments, a total of 3 ϫ 10 4 cells were plated in 96-well plates, incubated overnight, and treated with variable concentrations of PLX4032 for 60 min or 3.3 M PLX4032. Treatments were performed in biological triplicates. Plates were read on an En-Vision plate reader (PerkinElmer Life Sciences).
Cell Culturing of HCT116 and Colo205 Cell Lines-Cell lines were grown in double SILAC SD Media (Thermo Fisher Scientific) containing 10% FBS, 164 M lysine, 95 M arginine, and 4.3 M proline (Silantes, Munich, Germany) with additional nutrients consistent with Bendall et al. (28). Cell plates were incubated at 37°C and 5% CO 2 in a Panasonic incubator. Cells were counted using a Leica microscope with 10ϫ 0.25 objectives.
Approximately 12.5 ϫ 10 6 cells were plated in 32 125 ϫ 25-mm Petri dishes. Incubation with the Raf inhibitor vemurafenib (PLX4032) (Selleck Chemicals Houston, TX) was performed by adding 1 ml of PLX4032 or DMSO (vehicle) (Sigma-Aldrich) diluted in SILAC DMEM, 10% FBS to reach a final concentration of 3.3 M in each Petri dish. Cells were harvested every 5 min during the 1st h of treatment with either PLX4032 (heavy label) or DMSO (light label). Cells were collected by adding 10 ml of Ϫ80°C pre-cooled ethanol and final scraping. Heavy and light cell pellets were combined in 10-ml Falcon tubes. The workflow is shown in Fig. 1C.
Digestion and Desalting of Cell Extracts-Heavy and light cell pellets from each of the 15 time points were combined in 10-ml Falcon tubes, separated by centrifugation at 700 rpm (Sorvall Legend RT centrifuge), and washed with 5 ml of PBS. Cell lysis was performed by sonication for 15 s with a Sonic dismembrator (Thermo Fisher Scientific, Rockford, IL) after adding 1 ml of lysis buffer (8 M urea, 50 mM Tris, pH 8, HALT phosphatase inhibitor from Thermo Fisher Scientific). Cells were maintained at 0°C to prevent biological activities. Protein concentration was measured with Pierce BCA assay from Thermo Fisher Scientific. Protein disulfide bonds were reduced by incubating lysates in 5 mM dithiothreitol (DTT, Sigma-Aldrich) for 30 min at 56°C while shaking at 1000 rpm. Alkylation of cysteine residues was achieved by incubating iodoacetamide (Sigma-Aldrich) at a concentration of 15 mM for 30 min at room temperature in the dark. Excess iodoacetamide was deactivated by adding DTT up to 5 mM for 15 min at room temperature.
Phosphopeptide Enrichment-Phosphopeptide enrichment was performed on 5-m titansphere bulk particles (Canadian Life Science, Peterborough, Ontario, Canada) according to recent protocols (29,30). This method provides the most efficient enrichment of phosphopeptides.
Off-line Strong Cation Exchange Chromatography-To achieve high reproducibility and parallel sample fractionation, we used strong cation fractionation on homemade SCX spin columns packed with 18 -22 mg of polysulfoethyl A (300 Å bulk material, Canada Life Science; Peterborough, Ontario, Canada). After equilibrating the SCX spin columns with each 100 l of SCX buffer A (10% ACN, 0.2% FA), buffer B (1000 mM NaCl in 10% ACN, 0.2% FA), and finally 200 l of SCX buffer A, peptides were resuspended in 100 l of buffer A, loaded on the SCX column, and eluted in six salt step fractions of 100 l each (0, 30, 50, 80, 120, and 500 mM NaCl) in SCX buffer A. In subsequent analyses, the last salt step fraction was not injected as they contained negligible amounts of peptides. Prior to LC-MS/MS analyses, all fractions were dried on a SpeedVac centrifuge at room temperature (Thermo Fisher Scientific) and resuspended in 4% FA. All centrifugation steps were performed at 4°C except when indicated.
LC-MS/MS Analysis-All LC-MS/MS analyses were performed on a Q-Exactive plus mass spectrometer using homemade capillary LC columns (18 cm long, 150-m inner diameter, 360 m outer diameter). Capillary LC columns were packed with C18 Jupiter 3-m particles (Phenomenex, Torrance, CA) at 1000 p.s.i. Precolumns of 5 mm were packed at 100 p.s.i. using the same material. LC separations were performed at a flow rate of 0.6 l/min using a linear gradient of 5-40% aqueous ACN (0.2% FA) in 150 min. MS spectra were acquired with a resolution of 70,000 using a lock mass (m/z: 445.120025) followed by up to 10 MS/MS data-dependent scans on the most intense ions using high energy dissociation. AGC target values for MS and MS/MS scans were set to 1 ϫ 10 6 (maximum fill time 100 ms) and 5e4 (maximum fill time 120 ms), respectively. The precursor isolation window was set to m/z 1.6 with a high collision dissociation normalized collision energy of 24. The dynamic exclusion window was set to 30s.
Data Processing and Analysis-Raw data analysis of SILAC experiments was performed using MaxQuant software 1.3.0.5 and the Andromeda search engine (31). The false discovery rate (FDR) for peptide, protein, and site identification was set to 1%, the minimum peptide length was set to 6, and the "peptide requantification" function was enabled. The option match between runs (3-min time tolerance) was enabled to correlate identification and quantitation results across different runs. Up to two missed cleavages per peptide were allowed. To adjust for any systematic quantification errors, SILAC ratios were normalized by time point. All additional parameters are reported in MaxQuant parameters.txt and experimentalDesign.txt provided in supplemental Table S4. The Uniprot human proteome database (January, 2014, containing 70,584 entries) was used for all searches. Protein groups were formed by MaxQuant, and all identified peptides were used to make pairwise comparison between proteins in the database (31,32). Proteins containing an equal or overlapping set of peptides were merged together to form a protein group and ranked according to the highest number of peptides.
All raw LC-MS/MS data and MaxQuant output files can be accessed from PeptideAtlas using the identifier PASS00897 (at the following URL: PASS00897:PM8586nrx@ftp.peptideatlas.org/). Max-Quant output files contain lists of all matched and unmatched spectra, all matched protein isoforms and peptides including scorings and more detailed information.
In addition to an FDR of 1% set for peptide, protein, and phosphosite identification levels, we used additional criteria to increase data quality. We selected only peptides for which abundance ratios (FC ϭ PLX4032/Control) were measured in at least 8 of 13 time points. Then we set a cutoff for maximum phosphosite localization confidence across experiments (time points) to 0.75 (33). Next, we distinguished dynamic from static phosphosites by calculating an FDR based on curve fitting. This was achieved by fitting each profile to a polynomial order ranging from 1 to 8 and comparing the resulting profiles with those obtained from the same data following five rounds of randomization (supplemental Fig. S2). The coefficient of determination (R 2 ) was calculated for each profile of a given order and compared against those from the randomized dataset. For each polynomial order, dynamic profiles were selected based on R 2 values below an FDR of 1%. Next, we used fuzzy c-means algorithm to cluster all dynamic profiles (34). Analysis and visualization were performed in the R environment, and clusters were obtained with the Mfuzz package (35). Optimal setting of the "fuzzifier" parameter was 1.235 as estimated with the "mestimate" function. To find the optimal number of clusters, we performed repeated soft clustering for cluster numbers ranging from 2 to 30 and calculated the minimum centroid distance (minimum distance between two cluster centers produced by c-means clustering) (36). Next, we determined the membership value of dynamic profiles and selected only those with values greater or equal to 0.9 for high membership profiles.
Gene Ontology (GO) enrichment analyses were performed using the database DAVID version 6.7 (37). Protein interaction networks were defined using the STRING database (38), and interacting proteins were visualized with Cytoscape version 2.8.3 (39). Kinase substrate interactions were predicted using NetworKIN and Kinome-Xplorer (40). A protein-protein interaction network was built in STRING version 9.1 for all proteins containing dynamic phosphosites (38). All interaction predictions were based on experimental evidence with a minimal confidence score of 0.9 (considered as a "highest confidence" filter in STRING). Structure rendering was done using PyMOL (Schrö dinger).
Protein Purification and in Vitro Kinase Assays-GST-TEV fusion of human MKL1 (sequence corresponding to accession 3TPQ_M), MKL2(46 -162) (NP_054767.3), and TEAD3(1-436) (BC091488.1) were cloned between BamHI and XhoI in a pGEX4T1 plasmid engineered with a TEV cleavage site. The resulting plasmids were transformed in chemocompetent Escherichia coli BL21pLys, and expression was induced overnight at 18°C with 1 mM isopropyl 1-thio-␤-D-galactopyranoside in 500 ml of terrific broth. Cells were harvested and lysed by one cycle of freeze-thawing on ice in 25 ml of buffer F (10 mM Tris-HCl, pH 7.5, 150 mM NaCl, 1 mM EDTA, 1 mM PMSF) supplemented with benzonase (Sigma-Aldrich). Lysates were then sonicated four times for 10 s with a sonicator probe set at 80% power and cleared by centrifugation at 4°C for 20 min at 40,000 ϫ g in a Sorvall SS-34 rotor. The supernatant was adjusted to 0.05% Tween 20 and incubated with 150 l of glutathione-Sepharose 4B (GE Healthcare, Little Chalfont, UK) at 4°C for 4 h. Glutathione-Sepharose column was washed with 20 bed volumes of buffer G (50 mM Tris-HCl, pH 7.5, 50 mM NaCl, 0.1 mM CaCl 2 , 2 mM 2-mercaptoethanol, 0.05% Tween 20) and eluted at room temperature in buffer G with 25 mM reduced glutathione (pH was adjusted to 7.5 before elution). Recombinant protein concentration was established using the Bradford assay, and protein purity was confirmed by SDS-PAGE followed by Coomassie staining. Proteins were stored at Ϫ80°C in buffer G supplemented with 20% glycerol.
Active ERK was obtained from Carna Bio (Kobe, Japan) (04-143), and purified rabbit muscle actin was obtained from Cytoskeleton (AKL99). Actin was resuspended at 10 mg/ml in buffer H (5 mM Tris-HCl, pH 8.0, 0.2 mM CaCl 2 , 0.2 mM ATP, 5% sucrose). G-actin was prepared by incubation with 50 M latrunculin B for 1 h at room temperature followed by sedimentation of F-actin in a Sorvall Discovery M150 SE at 100,000 ϫ g for 1 h at 4°C.
In vitro kinase assays were conducted in kinase A buffer (50 mM HEPES, pH 7.5, 10 mM MgCl 2 , 1 mM EGTA, 0.01% Brij-35, 2 mM DTT, 1 mM PMSF, 1ϫ leupeptin and pepstatin) with 30 ng of active ERK and 2 g of GST fusion protein. When required, the reactions were combined with TEV cleavage of GST fusions by including 3 units of acTEV (Invitrogen) in the kinase reactions. Kinase reactions were initiated by addition of ATP (50 M final) ϩ 5 Ci of [␥-32 P]ATP (PerkinElmer Life Sciences) and were incubated at 37°C for 30 min. Reactions were stopped at 95°C in the presence of Laemmli sample buffer. Kinase reactions were loaded on SDS-polyacrylamide gels; proteins were stained with Coomassie G-250, and gels were finally dried and exposed to X-ray films.
For actin titration experiments, recombinant MKL1 was incubated overnight with different molar ratios of G-actin. The protein mixes were then subjected to kinase reaction in the presence of ERK2 and ATP as described above. All experiments were performed in biological triplicates, and only the best Western blottings are depicted in the corresponding figures.

RESULTS
Dynamic Changes in the Phosphoproteome Following Vemurafenib Treatment-For this study, we chose two colon cancer cell lines that display elevated ERK signaling due to different genetic lesions. Colo205 cells harbor a BRAF V600E mutation but are otherwise WT for RAS, whereas HCT116 cells contain a KRAS G13D mutation but are WT for RAF (41). These cell lines are widely used in the cancer research community and respond differently to RAF inhibitors. Indeed, phospho-ERK1/2 (pERK) is inhibited in Colo205 cells following incubation with the RAF inhibitor, although it is induced in HCT116 cells (20,42). We profiled the change in ERK1/2 activation for different concentrations of the RAF inhibitor vemurafenib, and we determined that a concentration of 3.3 M induced an 8-fold increase of pERK in HCT116 cells (Fig.  1A). In contrast, the same concentration resulted in full pERK inhibition in Colo205 cells (Fig. 1A). Time course measurements showed that the pERK response to vemurafenib treatment was maximal after 20 min in both cell lines (Fig. 1B).
To evaluate the downstream effects of this response, we monitored the global change in the phosphoproteome of both cell lines over a period of 60 min post-stimulation with a temporal resolution of 5 min (Fig. 1C). We used SILAC (43) where control (DMSO vehicle) and stimulated cells (3.3 M vemurafenib) were cultured in dialyzed FBS for more than 10 doublings with light and heavy Lys and Arg amino acids (see under "Material and Methods" for additional details). Cells from each time point were harvested, rinsed with PBS, snapfrozen in precooled ethanol at Ϫ80°C, and then combined in equal amounts. Prior to MS analyses, cell suspensions were washed with PBS to remove ethanol and lysed in urea, and proteins were reduced, alkylated, and subsequently digested with trypsin. Phosphopeptides were enriched using TiO 2 microcolumns and fractionated off line into six SCX salt fractions prior to LC-MS/MS analysis. Phosphopeptide identification and quantification were performed using MaxQuant and processed in R to cluster phosphopeptide profiles and determine regulated phosphorylation sites.
These analyses enabled the identification of 28,711 and 20,931 phosphopeptides (localization confidence Ն75%) for HCT116 and Colo205 cells, respectively. In total, we identified 37,910 unique phosphopeptides, of which 10,300 were quantified across more than eight time points (supplemental Fig. S1). To profile the global changes in the phosphoproteome, we determined the FC in phosphopeptide abundance for the 1st h after incubation with vemurafenib ( Fig. 2A).
To define regulated phosphosites in the dataset, we first selected phosphopeptides quantified in more than eight time points and performed fitting of all kinetic profiles with a polynomial model. Next, we distinguished dynamic from static phosphosites by calculating a false discovery rate based on curve fitting. This was achieved by fitting each profile to a polynomial order ranging from 1 to 8 and comparing the resulting profiles with those obtained from the same data following five rounds of randomization (supplemental Fig. S2). The coefficient of determination (R 2 ) was calculated for each profile of a given order and compared against those from the randomized dataset. For each polynomial order, dynamic profiles were selected based on R 2 values below an FDR of 1% and prevented over-fitting profiles to high polynomial order. Accordingly, these analyses identified 660 dynamic phosphopeptide profiles out of 10,300 profiles for polynomial order Յ4 (supplemental Table S1). These dynamic profiles formed four distinct clusters according to the Mfuzz algorithm (35), and were either up-regulated or down-regulated and showed an early or a late response (Fig. 2C). From the 660 dynamic profiles, we selected only those that have membership values equal or greater than 0.9, which resulted in 220 "high membership" profiles (supplemental Table S2).
Dynamic data show progressive broadening of FC values observed for both cell lines with a gradual shift in the median distribution toward a decrease or an increase in phosphorylation for Colo205 and HCT116, respectively. The width of this distribution was calculated from the interquartile range (IQR) for the log 2 (FC) distribution of phosphopeptides at each time point and plotted over the entire cell stimulation period (Fig.  2B). The resulting curve showed a progressive and linear increase in the phosphopeptides fold-change IQR with time, a trend that was markedly different from non-phosphorylated peptides that remained largely unaffected over the same time period (Fig. 2B). Also, changes in protein phosphorylation were more pronounced for HCT116 compared with Colo205 as reflected from the variation of IQR with time. Interestingly, we observed a steady increase in IQR values during the 1st h, although phospho-ERK1/2 activation/inhibition levels reached a maximum after 20 min post-stimulation (Fig. 1B).

RAS-ERK Pathway Modulation Accounts for Most Dynamic Phosphoproteome Changes Induced by Vemurafenib-Pro-
files within the high membership dataset were used to contrast the response of the two cell lines and to predict kinase activities following vemurafenib inhibition. The grouping of profiles from both cell types identified five pairs of profiles with distinct patterns (Fig. 3A). The first group is represented by opposite temporal changes in the two cell lines and cor-respond to 33 down-regulated phosphosites identified in Colo205, and 32 were up-regulated in HCT116 (Fig. 3A). The second and most populated group comprised 124 phosphopeptides up-regulated in HCT116 and were either not identified (92) or showed a static behavior in Colo205 (32) (Fig. 3A). The third group corresponds to a more limited number of phosphopeptides that were dephosphorylated in Colo205 but were either not identified (24) or remained static in HCT116 (2) (Fig. 3A). Groups 4 and 5 are discussed below.
The first three groups shown in Fig. 3A represent 83% of the high membership sites and are directly or indirectly associated with RAS-ERK pathway activation (Fig. 3A, ERK-like profiles). Consistent with this proposal, phosphopeptides comprising MEK2 (MAP2K2) and ERK1 (MAPK3) activation sites were clustered in groups 1 and 2, respectively (Fig. 3, A  and B). A large proportion of dynamic phosphorylation sites thus behaved in accordance to the RAS-ERK pathway response, whereby up-regulated phosphopeptides were largely dominated by profiles identified in HCT116 cells (156/165), whereas down-regulated phosphopeptides were mostly found in Colo205 (59/88). Also, the majority of the 220 dynamic phosphorylation sites of our high membership dataset comprised proline-directed phosphorylation motifs (148/220), a hallmark of ERK-mediated events (6). However, we noted that several sites of the high membership dataset were not fulfilling the ERK consensus motif and could represent secondary or tertiary events occurring within the highly inter-connected kinase-substrate network downstream of ERK. To predict which kinase-substrate relationships were responsible for dynamic phosphorylation events in our dataset, we used the bioinformatic tool KinomeXplorer (40), a platform based on an improved NetworKIN algorithm that models kinase signaling networks by combining linear motif analysis with protein-protein interactions. The corresponding analysis is shown in Fig. 3C where kinase groups are clustered according to KinomeXplorer. This analysis confirmed that a substantial enrichment of dynamic phosphorylation sites are predicted to be ERK-mediated (Fig. 3C). In addition, putative substrates of JNK and p38 comprised a proline-directed consensus motif and formed a subset of the dynamic sites (Fig. 3C). This analysis also indicated that our high membership dataset contained putative PKB, GSK3, PAK (p21-activated kinase), CK1, and p70 S6K substrates. Interestingly, two AGC kinases (PKB and p70 S6K ) share a RXRXX(S/T) consensus motif with known ERK downstream effectors of the RSK family (p90 RSK or RPS6KA) (44). Concurrently, the RSK paralogs RPS6KA1, -3, and -4 were dynamically modulated on ERK-dependent sites in our dataset (supplemental Table S1), and therefore, several dynamic phosphorylations on RXRXX(S/T) sites are likely to be secondary events associated with RSK proteins. Supporting this idea, several known RSK-modulated phosphosites (PDCD4_ S457, BAD_S75, and TBC1D4_S588) were identified in our dataset (supplemental Tables S1 and S2) (45)(46)(47).

Light Heavy
A small subset of 28 sites in HCT116 and 9 sites in Colo205 from our high membership dataset did not correlate with the pERK response and could potentially represent off-target events (Fig. 3A, groups 4 and 5). A closer examination of these data revealed that several sites are located on peptides that contained multiply phosphorylated forms. The modulation of these phosphopeptides with a kinetic profile inconsistent with the ERK activation/inhibition status is explained in part by the processivity of phosphotransfer whereby we observed an increase in the abundance of the multiphosphorylated peptides, although the singly phosphorylated form decreased in abundance. The overall level of phosphorylation of a given peptide can therefore increase despite the fact that the monophos- phorylated form disappeared over time. A case in point is ETV3 whose phosphorylation on Ser-245 rapidly decreased, whereas the doubly phosphorylated form (Ser-245/ Ser-250) increased during the same time period in HCT116 (Fig. 3B). It is noteworthy that 23 phosphopeptides from group 4 are singly phosphorylated on a proline-directed motif, of which 13 contained a second putative proline-directed phosphorylation site (supplemental Table S3). Altogether, these observations support the notion that most activity associated with vemurafenib treatment corresponds to on-target effects with a limited number of secondary and tertiary signaling events associated with downstream ERK-modulated substrates.  These observations are consistent with a recent phosphoproteomics study comparing RAF and MEK inhibitor activity (26).

A)
Time-resolved Phosphoproteome Profiling Identifies New ERK Substrates-Our analysis established that dynamic phosphorylation sites were mostly associated with the regulation of the RAS-ERK pathway. We next wanted to exploit the kinetic information in our dataset to identify novel high affinity downstream substrates of ERK. We posit that the clustering of phosphopeptides as "Early" or "Late" could reflect ERKsubstrate affinities. To perform this analysis, we grouped all high membership phosphorylation sites from clusters "Early Down" and "Early Up" in both cancer cell lines and searched for ERK1/2 physical interactions with cognate proteins using the STRING databases (48). This network was complemented with information on known and predicted ERK substrates derived from the BioGRID database (49) and the NetworKIN algorithm (50), respectively. This enabled the identification of 22 uncharacterized interactions between ERK and putative Early substrates (Fig. 4, red and yellow nodes). The similarity of these putative ERK-mediated kinetic profiles was displayed as the average change in phosphorylation over time (Fig. 4, top right).
We complemented this interaction-based analysis with predictions of docking motifs to determine whether ERK selectivity and processivity was also influenced by sequence motifs present in the corresponding substrates. We used Scan-Prosite (51) to identify KIM motif (D-motif; variable consensus) (6) and DEF motif (FXFP) (52) within 50 amino acid distance of the high membership ERK phosphorylation sites (Fig. 4, blue dots, and supplemental Fig. 3A). We observed an enrichment of KIM and DEF motif consensus sequences in our dynamic set compared with static phosphoproteins, although only a small subset of 13 sequences included these docking sites. Furthermore, six high membership HCT116 phosphorylation sites that contained a docking motif were grouped in the Early cluster (supplemental Fig. 3B). In contrast, KIM-or DEFcontaining proteins were distributed in both Early and Late clusters in Colo205 cells (supplemental Fig. 3B). This result is not entirely unexpected because dephosphorylation of ERK substrates in Colo205 does not solely rely on ERK-substrate affinity but also on phosphatase-substrate interactions. Interestingly, the average kinetic profiles of Early phosphosites with KIM or DEF motifs showed a more rapid change in phosphorylation compared with those that did not contain these motifs (supplemental Fig. 3C). Examples of phosphosites for substrates with and without KIM or DEF motifs are shown in supplemental Fig. 3, D and E. It is noteworthy that bioinformatic searches for sequence determinants other than KIM and DEF motifs failed to identify relevant linear motifs in our datasets that could play a role in ERK-substrate interactions (data not shown).

Time-resolved Phosphoproteome Analysis Identifies Novel Biological Functions Connected with RAS-ERK Signaling-To
decipher which cellular functions were regulated upon RAS-ERK pathway modulation by vemurafenib, we performed standard GO enrichment analysis of our high membership and dynamic datasets using DAVID (37,53). These analyses indicated the enrichment of genes associated with "RAS/ RTK signaling," "regulation of Rho signal transduction," "kinases/phosphatases," and "actin cytoskeleton organization" (Table I Not Enriched). The dynamic gene set yielded similar GO term enrichment but also included the functional groups related to the regulation of "cell death/apoptosis," "mRNA transport/nuclear pore," and "RNA splicing" (Table I). Although several of these GO terms were previously shown to connect with the RAS-ERK pathway, such as mRNA transport/nuclear pore (54), our analysis identified other less usual connections between ERK and "actin cytoskeleton organization" and "regulation of Rho signal transduction." We used GO terms enriched in our dataset (Table I) combined with manual curation of gene lists to derive a systems level view of these novel RAS-ERK pathway connections (Fig. 5).
These analyses revealed that several phosphosites occurred upstream of RAF and may represent negative feedback loops. Among those, we identified the time-dependent modulation of Ser-642 in CRAF (RAF1) C-terminal tail, a target site of activated ERK involved in desensitizing CRAF to additional stimuli (55). Upstream factors of canonical RTK-RAS signaling (RTKs ERBB2 and EGF receptor and the RAS-GEF SOS1) as well as the insulin receptor (InsR) signaling pathway (IRS1/2, GIGYF1/2) were also phosphorylated at proline-directed sites. Similarly, scaffolding proteins of the RAS-ERK cascade (GRB10, PXN, and CNKSR1) were dynamically phosphorylated in response to vemurafenib treatment (Fig. 5). Interestingly, several of these putative feedback sites were modulated early after vemurafenib treatment, suggesting that phosphorylation of these sites is concurrent with the phosphorylation of downstream ERK effectors.
The phosphorylation of nuclear pore complex (NPC) subunits by ERK was reported in an earlier phosphoproteomic study (54). Here, we identified the dynamic phosphorylation of NUP50, NUP107, NUP153, NUP188, TPR, and RANBP2 (Figs. 5 and 6 and supplemental Table S1). These proteins participate in several sub-structures of the nuclear pore and have distinct functions during translocation of biomolecules STRING was used to obtain high confidence ERK1/2 interactors among proteins with dynamic phosphosites. This was complemented with known ERK1/2-regulated sites reported in BioGRID. Dynamic phosphosites were found in known interacting proteins (red dots, StringDB), in predicted direct interactors (yellow dots, Networkin algorithm), or in validated ERK substrates (green dots, BioGRID interaction database). Nodes represented in this high confidence ERK substrate network display similar dynamic behaviors upon vemurafenib treatment (top right, shaded area across averaged profile represents 95% confidence interval). across the pore. ERK-dependent phosphorylation of NPC proteins was shown to affect the individual properties of the subunits. For instance, NUP50 phosphorylation on ERK sites located within FG repeats regulates its interaction with importin-␤ and transportin (54). Similarly, mutation of ERK sites in TPR resulted in cell cycle defects associated with the failure to recruit MAD1/2 to the NPC (59).
ERK activation is necessary for its migration into the nucleus where it phosphorylates specific regulatory proteins and transcription factors to initiate transcription of immediateearly genes (60 -62). Accordingly, quantitative phosphoproteomics identified several known transcription factors phosphorylated by ERK such as ETV3_S245/S250, ERF_S531, JUND_S90/S255, and STAT3_S727 (63, 64) as well as putative substrates such as FOXP4_S554, TP53BP1_S366, and TEAD3_S148 (Figs. 5 and 6). TEAD3 participates in the Hippo tumor-suppressor pathway, which controls organ size and cell death. This pathway culminates in the phosphorylation and inactivation of the transcriptional co-activator YAP1, which itself relies on DNA-binding adaptors of the TEAD family to mediate a complex transcriptional response (65). YAP1_S61 is a bona fide substrate of the Hippo pathway kinase LATS1 (66) and was found up-regulated only in Colo205 cells. In contrast, TEAD3_S148, a proline-directed site, showed an increase in phosphorylation upon vemurafenib treatment only in HCT116 cells. Mammalian TEAD proteins (TEAD1-4) share a common domain architecture consisting of a TEA DNA-binding domain and a YAP1interacting region (67). TEAD3_S148 is located in the linker region between these two conserved portions of the protein.
To assess direct phosphorylation of this site by ERK, we performed in vitro kinase assays. Although WT TEAD3 was phosphorylated in the presence of active ERK2, 32 P incorporation was completely abrogated with TEAD3 S148A variant (Fig. 6). This establishes a direct connection between RAS-ERK signaling and the Hippo pathway transcription factor TEAD3.
Among the transcription factors dynamically modulated upon vemurafenib treatment, our attention was also attracted by the presence of megakaryoblastic leukemia 1 (MKL1) and 2 (MKL2). These were phosphorylated on Ser-33 and Ser-77, respectively, two orthologous phosphosites that comprise an upstream KIM motif and that were not previously characterized (Figs. 4 and 5). Although these sites were not identified as ERK substrates based on STRING interactions, a previous report suggested the direct phosphorylation of MKL1_S454 by ERK as part of a negative feedback mechanism (68). Interestingly, our phosphoproteomic analyses revealed the rapid phosphorylation of MKL1_S33 and MKL2_S77 upon vemurafenib treatment, whereas MKL1_ S454 remained unchanged ( Fig. 6 and supplemental Fig. S4). Further experimental validation was obtained using in vitro kinase assays. GST-MKL fusion substrates were incubated with active ERK in the presence of radiolabeled ATP. In vitro reactions confirmed that active ERK specifically phosphorylated MKL1_S33 and MKL2_S77 because 32 P incorporation was abrogated in S33A and S77A mutant proteins (Fig. 6).

MKL1 Ser-33 Defines a Novel Actin-dependent Mechanism That Modulates Kinase-Substrate
Interactions-MKL1/2 were grouped under the GO function "actin cytoskeleton organization," a term being highly enriched in our dataset. These proteins act as actin-dependent transcriptional co-activators of the serum-response factor (SRF), a known transcriptional regulatory hub acting downstream of the RAS-ERK pathway (69). In resting cells, MKL1 transit between the nucleus and the cytosol, and its localization was reported to depend on the availability and binding of globular actin (G-actin) to RPEL motifs. In serum-starved conditions, MKL1 is bound to G-actin and sequestered in the cytoplasm but is rapidly translocated to the nucleus in conditions that favor the polymerization of actin and the formation of stress fibers (70 -73). Structural analyses of MKL1-actin interaction revealed that RPEL motifs and the two intervening spacers can each bind a G-actin monomer, forming a compact assembly of five G-actins with one MKL1 RPEL repeat (Fig. 7A) (74 -76). Close inspection of MKL1 crystal structures revealed that KIM residues overlap with the G-actin-binding site found in RPEL1 (Fig. 7, A and B, green) (71,73). Furthermore, MKL1 Ser-33 is tightly packed against the second G-actin mole-

TABLE I Gene ontology terms (biological function) enriched among phosphorylated proteins that show time-dependent modulation after vemurafenib treatment
Comparison of the MFuzz high membership dataset (clustering membership Ն0.9, total of 220 phosphorylation profiles) and dynamic dataset (clustering membership Ն0.0, total of 660 phosphorylation profiles) is shown. cule that interacts with the linker connecting RPEL1 and RPEL2 (Fig. 7A, yellow). Based on these observations, we posit that the binding of G-actin to RPEL1 domain of MKL1 (green) might control the access of ERK to the conserved KIM motif, thereby affecting the phosphorylation of MKL1 Ser-33 (Fig. 7C). Also, MKL1 Ser-33 must be released from the second G-actin subunit (Fig. 7C, yellow) to facilitate its phosphorylation by ERK. To validate this proposal, we performed in vitro kinase assays on recombinant MKL1 harboring KIM motif mutants R5A/R6A and L9A/L11A or RPEL1 mutant R16A/R17A (Fig. 7D). These analyses confirmed that MKL1 phosphorylation by ERK was abrogated for the R5A/R6A and L9A/L11A MKL1 mutants, suggesting that the KIM motif is necessary for docking onto ERK to facilitate MKL1 Ser-33 phosphorylation. In contrast, the R16A/R17A double mutant known to affect MKL1 binding to G-actin (74) did not affect ERK-dependent phosphorylation of Ser-33 (Fig. 7D).
To determine whether G-actin binding to MKL1 could affect its phosphorylation by ERK, we preincubated WT FIG. 5. Depiction of RAS-ERK signaling pathway and primary and secondary connections to substrates identified by time-resolved phosphoproteome profiling. Primary and secondary targets dynamically modulated upon vemurafenib treatment. Main functional groups depicted are "RPS6 kinases," "RTK signaling," actin cytoskeleton organization, nuclear pore, and "transcription." Key players of actin cytoskeleton regulation, including Rho GEF's and GAP's, are emphasized. The Hippo pathway target, YAP1, acts together with TEAD factors (here TEAD3) to control transcription at promoters containing M-CAT elements. TCF and MKL1/2 assemble with SRF on promoters comprising serum response elements (SREs). The MKL1/2 regulators LMO7 and FLNA/B included putative RSK (RPS6A) phosphosites that were dynamically modulated after vemurafenib treatment.
MKL1 overnight with various concentrations of G-actin and performed in vitro ERK assays (Fig. 7E). These experiments revealed that phosphorylation of MKL1 was only observed at low G-actin concentration, suggesting that excess actin shielded MKL1 Ser-33 and prevented its phosphorylation by ERK (Fig. 7E, left panel). To rule out a possible nonspecific interaction due to protein crowding in the kinase assay, we conducted the same in vitro kinase assay using the MKL1 R16A/R17A double mutant (Fig. 7E, right panel). In contrast to WT MKL1, the R16A/R17A double mutant was fully phosphorylated irrespective of G-actin concentration. These experiments support the notion that actin binding to MKL1 prevents Ser-33 phosphorylation by ERK and suggest an interplay between actin dynamics and MKL phosphorylation. In the course of preparing this paper, consistent findings were published by Treisman and co-workers (77), who validated Ser-33 phosphorylation by ERK as well as the importance of the upstream KIM motif and the role of globular actin binding in blocking phosphorylation at this site. However, we note that residue numbering differed between this publication (Ser-98) and ours (Ser-33) because different mRNA isoforms were used as references (Ser-33 numbering corresponds to MKL1 isoform X1 (XP_ 016884377)).

DISCUSSION
Mass spectrometry combined with metabolic labeling (e.g. SILAC) enabled the quantification of phosphoproteins on a systems-wide scale. The temporal profiling of protein phosphorylation facilitated the identification of RAS-ERK pathway substrates and secondary events in two colon cancer cell lines that respond differently to the RAF inhibitor vemurafenib. By taking advantage of the opposite response to RAF inhibitors, our experimental strategy provided a quantitative view of ERK signaling dynamics without the pleiotropic effects of growth factors, serum, or phorbol esters that can give rise to pathway cross-talks. Our quantitative phosphoproteomic analyses revealed that vemurafenib has low off-target activity on kinases other than RAF or on phosphatases. This observation is consistent with a recent study comparing the proteome-wide effects of BRAF and MEK inhibitors when used alone or in combination (26).
The use of targeted RAF stimulation in combination with time-resolved phosphoproteomics facilitated the deconvolution of early and late phosphorylation events and associated downstream and upstream targets of ERK. Our analysis suggests that predicted direct ERK substrates (i.e. containing a  (76). B, schematic representation of a model whereby actin binding to MKL1 prevents phosphorylation of MKL1 Ser-33 by ERK1/2. C, sequence alignment of MKL1 showing the conservation of KIM motif and Ser-33 across species. Residue Ser-33 is located between the first two RPEL repeats required for monomeric actin binding. The putative kinase-interaction motif (KIM) for ERK as well as the actin-binding region of the first RPEL motif are also highlighted. D, in vitro kinase assays of wild type and alanine mutant GST-TEV-MKL1 with recombinant active ERK2. Mutation of residues Arg-5 and Lys-6 or Leu-9 and Leu-11 within the KIM motif is sufficient to prevent phosphorylation of Ser-33 by ERK2. E, in vitro phosphorylation assay of wild type and alanine mutant MKL1 in presence of recombinant active ERK2 and variable amounts of actin. Phosphorylation of MKL1_S33 by ERK2 is impeded at elevated actin concentrations (8:1 and 16:1 molar ratios), whereas the phosphorylation of R16A/R17A mutant MKL1 is resistant to actin titration. proline-directed site) are in general rapidly phosphorylated upon RAF induction, as evidenced by the early phosphorylation of p90 RSK proteins RPS6KA1 and RPS6KA3. A rapid induction was also observed for other well known substrates such as nuclear pore complex subunits and transcription factors (e.g. ETV3, JUND, and ERF). Rapid kinetics was also observed for uncharacterized ERK downstream targets like several Rho GEFs and Rho GAPs. Nonetheless, bona fide ERK target sites were also associated with late events (e.g. STAT3, NUP188, and RPS6KA4). Therefore, substrate phosphorylation kinetics alone cannot distinguish direct from indirect ERK targets. Kinetic differences might, however, be explained by distinctive molecular features such as linear motifs (DEF and KIM motifs), three-dimensional structure elements, quaternary complex formations, and cellular co-localization that could account for changes in kinase-substrate affinities.
In addition to ERK targets, our study provided insights into the dynamic changes in phosphorylation of secondary events such as the phosphorylation of several RSK substrates. Interestingly, several known (TBC1D4, BAD, and PDCD4) and predicted (FLNA and FLNB) RSK target sites were identified as late events in our dynamic phosphoproteomic dataset (supplemental Tables S1 and S2). This observation is consistent with the view that RSK-dependent signaling events are subsequent to and thus downstream of ERK signaling. Although this will require more data and analysis, it is likely that time-resolved phosphoproteomics will allow epistatic analysis of cell signaling networks.
Dynamic phosphoproteomics identified several known (e.g. STAT3, ETV3, ERF, and JUND) and putative ERK substrates (e.g. TEAD3, TAF6, TCF12, FOXK2, and FOXP4) involved in transcriptional regulation. In particular, we demonstrated for the first time that the transcription factor TEAD3 is directly phosphorylated by ERK at Ser-148 in the linker region between TEA DNA-binding domain and YAP1-interacting region (67). Our results are thus consistent with the presence of cross-talk between the RAS-ERK and the Hippo pathway that would implicate TEAD transcription factors. The phosphorylation of TEAD3 at Ser-148 could modulate its interaction to DNA and/or to YAP1 either alone or through the recruitment of other binding partners. YAP1 is a downstream regulatory target in the Hippo signaling pathway, and its binding to TEAD transcription factors is required to stimulate the expression of genes involved in tissue homeostasis (78). Interestingly, a recent genetic screen in BRAF mutant tumor cells identified YAP1 as a survival input mediating RAF and MEK inhibitor resistance in melanoma (79). This, together with other evidence (80, 81), suggests a cross-talk between the Hippo and RAS-ERK pathways. The role of TEAD3 phosphorylation on Ser-148 in this cross-talk would deserve to be investigated.
Our study also identified the RPEL motif of MKL transcription factors as putative ERK substrates. Using in vitro kinase assays we demonstrated that ERK directly phosphorylates MKL1_S33 and MKL2_S77 in the linker region between RPEL1 and RPEL2 domains. The phosphorylation of this orthologous site in both proteins relies on a canonical ERK KIM motif that was only recently identified (77). Our analysis showed that this motif facilitates the binding of ERK to MKL1 and enhances its phosphorylation. We also confirmed that accessibility to this KIM motif and to Ser-33 is modulated by G-actin binding to the MKL1 RPEL repeats, consistent with results described in Panayitou et al. (77). On the basis of these results, we proposed a novel mode of kinase-substrate regulation whereby G-actin binding can mask the docking and phosphorylation of MKL1/2 by ERK. Alternatively, unbound MKL1/2 can be more readily phosphorylated by ERK when actin dynamics is shifted to the formation of filaments. It will be interesting to determine whether actin-dependent regulation of phosphosites applies more widely to other kinasesubstrate pairs. Conversely, it was shown by gel filtration that Ser-33 phosphorylation has an impact on the stoichiometry of actin binding by MKL1 RPEL motifs (Ser-98 in reference 77). Indeed, the WT protein binds to three molecules of actin, whereas the phosphomimetic mutant binds to two actin molecules (77). This result implicated ERK-dependent phosphorylation in the regulation of MKL1 nucleocytoplasmic shuttling by globular actin, whereby a reduction in actin binding upon Ser-33 phosphorylation promotes nuclear import (77).
Interestingly, our analysis identified several Rho regulators presumably phosphorylated by ERK upon vemurafenib treatment (ARHGAPs and ARHGEFs). It is unclear whether these phosphorylation events in ARHGAPs and ARHGEFs affect their activity and downstream Rho activation status. Nevertheless, it is possible that these factors indirectly impinge on MKL1/2 activity by modulating the balance between globular and filamentous actin. Our analysis additionally identified LMO7 and FLNA as potential RSK substrates. These are two recently reported regulators of MKL1, further reinforcing the role of ERK in MKL1/2-dependent transcriptional control. LMO7 was shown to relieve actin-mediated repression of MKL, although a recent study indicated that FLNA directly interacts with MKL1 (82,83).
It is noteworthy that unlike another study conducted by Muehlich et al. (68), our data do not provide evidence for feedback phosphorylation sites Thr-450 and Ser-454 on MKL1/2 upon specific RAF-MEK-ERK activation. This might be explained either by the different stimulation periods (e.g. 0 -60 min of vemurafenib treatment versus up to 3 h by Muehlich et al. (68)) or by the different stimulus used to trigger activation of the pathway. For instance, serum stimulation used by Muehlich et al. (68) triggers a pleiotropic response that includes RAS and Rho activation, whereas vemurafenib acts specifically downstream of RAS to selectively turn on the RAF-MEK-ERK signaling pathway. Alternatively, feedback phosphorylation of MKL1 may require elevated ERK activation, which might not be achieved when paradoxical activation of the pathway is induced by vemurafenib. This novel link between ERK signaling and MKL1/2 regulation suggests that at least two ERK outputs converge on cofactors of SRF (Fig. 5). The first one targets Glu-26 transformation-specific transcription factor family proteins (e.g. ETV3 and ERF, identified in this study), together known as ternary complex factors (TCFs), that associate with SRFbound promoters and are well known ERK substrates (84). The second one consists of MKL1/2 transcriptional co-activators. We thus suggest that TCFs might act in conjunction with an MKL1/2 axis downstream of ERK to modulate the transcription of SRF-bound elements (serum-response elements; CArG) in promoters (Fig. 5). This possible cross-talk between RAS-ERK and Rho-actin-MKL-SRF signaling is interesting on at least two accounts. First, RAF mutations (BRAF or CRAF) are found in Noonan syndrome and cardiofacio-cutaneous syndrome that belong to a group of syndromes called RASopathies but also in childhood-onset dilated cardiomyopathy. The clinical presentation of patients afflicted with these syndromes includes cardiac defects often associated with heart hypertrophy, which is also recapitulated in mouse models of Noonan syndrome (85,86). Myocardin family transcription factors, including myocardin itself (MYOCD), MKL1, and MKL2 are known to regulate an SRFdependent transcriptional program involved in heart development (87)(88)(89)(90)(91). It is thus possible that this RAS-ERK-MKL1/2 signaling axis contributes to some of the clinical manifestations observed in RASopathies. Second, several lines of evidence tie Rho signaling (92), MKL1, and SRF-dependent transcription with epithelial-mesenchymal transition and metastasis (93). The Rho-MKL1-SRF axis may therefore crosstalk with RAS-ERK signaling to modulate cancer cells invasiveness.
We anticipate that the global profiling of phosphorylation dynamics in response to cell stimuli will provide the foundation for the deconvolution of the cellular phosphoregulatory network. Phosphokinetic experiments such as those described here will also provide an important resource to further rationalize direct/indirect interactions and off-target effects between substrates and their putative kinases.
Data repository: The raw data arising from this project has been uploaded into the PeptideAtlas repository, with the Identifier: http://www.peptideatlas.org/PASS/PASS00897, project id: PASS00897.